From 584b500df42b78a0485a1c012345e8dbc3a99e60 Mon Sep 17 00:00:00 2001 From: Nick Newson Date: Sun, 13 Sep 2026 15:29:41 +0100 Subject: [PATCH 1/5] One comparison authority, and NaN stops being equal to everything MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every approximate comparison in the engine — VecBase, Mat3, Mat4, Quaternion — was the same four lines: const float diff = a - b; if (diff > eps || diff < -eps) { return false; } Both comparisons are false when diff is NaN, so a NaN compared approximately equal to anything, itself included. Two equal infinities subtract to NaN and passed the same way, which reads as correct until you notice that opposite infinities do too. A test asserting a transform had stayed finite would pass on a transform that was entirely NaN: precisely the failures tests exist to catch, hidden by the thing catching them. math/scalar.hpp is now the only place that decides whether two numbers are close. It asks `a == b` first, so identical values are equal whatever they are — including equal infinities, which the old form got right by accident rather than by rule. Any remaining non-finite operand is unequal. Finite values compare against the larger of an absolute and a relative tolerance, in double intermediates: `a - b` in float overflows for values near FLT_MAX of opposite sign, so the difference itself could be non-finite for two perfectly finite inputs and the comparison would answer about a number neither caller mentioned. THREE OVERLOADS, because a tolerance a caller spells out must mean what it says. No tolerance takes both defaults; one tolerance is ABSOLUTE ONLY, exactly the historical comparison; two are both terms, stated. A single function with a defaulted relative term would have made approxEqual(rhs, 1e-9f) admit a 1e-7 difference at magnitude 1.0 — an explicit tolerance silently overridden by an implicit one. The existing strictness tests are unchanged, which is the evidence that the overloads did their job. TOLERANCES ARE VALIDATED, not reinterpreted. A negative, NaN or infinite tolerance returns false, and the check runs before the `a == b` shortcut so equal operands cannot hide it. Treating a negative as "a term that cannot be satisfied" is mathematically defensible and diagnostically useless: it converts a broken configuration into a different valid policy and lets the run continue looking healthy. float_relative_epsilon is its own constant. float_epsilon was already serving as absolute tolerance and degeneracy threshold; a third job would have made it a number nobody could change safely. --- CMakeLists.txt | 1 + include/fire_engine/math/mat3.hpp | 25 +++- include/fire_engine/math/mat4.hpp | 25 +++- include/fire_engine/math/quaternion.hpp | 34 ++++-- include/fire_engine/math/scalar.hpp | 118 +++++++++++++++++++ include/fire_engine/math/vec_base.hpp | 25 +++- tests/math/test_scalar.cpp | 150 ++++++++++++++++++++++++ 7 files changed, 361 insertions(+), 17 deletions(-) create mode 100644 include/fire_engine/math/scalar.hpp create mode 100644 tests/math/test_scalar.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 855b4139..be8a6f47 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -410,6 +410,7 @@ add_executable(test_fire_engine tests/physics/test_physics_handle.cpp tests/physics/test_physics_determinism.cpp tests/physics/test_demos.cpp + tests/math/test_scalar.cpp tests/math/test_mat3.cpp tests/math/test_singular_value.cpp tests/math/test_mat4.cpp diff --git a/include/fire_engine/math/mat3.hpp b/include/fire_engine/math/mat3.hpp index f38ec2ad..cec2b39c 100644 --- a/include/fire_engine/math/mat3.hpp +++ b/include/fire_engine/math/mat3.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include namespace fire_engine @@ -227,13 +228,18 @@ class Mat3 return true; } + // Approximate equality, component by component, through the ONE scalar authority + // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit + // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart + // than 1e-9, exactly as it always did), and both arguments mean both terms. An invalid + // tolerance — negative, NaN or infinite — makes the comparison FALSE rather than being + // reinterpreted. NaNs compare unequal now, which is the defect this replaced. [[nodiscard]] - constexpr bool approxEqual(const Mat3& rhs, float eps = float_epsilon) const noexcept + constexpr bool approxEqual(const Mat3& rhs, float eps, float relativeEps) const noexcept { for (int i = 0; i < 9; ++i) { - const float d = m_[i] - rhs.m_[i]; - if (d > eps || d < -eps) + if (!almostEqual(m_[i], rhs.m_[i], eps, relativeEps)) { return false; } @@ -241,6 +247,19 @@ class Mat3 return true; } + [[nodiscard]] + constexpr bool approxEqual(const Mat3& rhs, float eps) const noexcept + { + // ABSOLUTE ONLY — a stated tolerance is the whole answer. + return approxEqual(rhs, eps, 0.0f); + } + + [[nodiscard]] + constexpr bool approxEqual(const Mat3& rhs) const noexcept + { + return approxEqual(rhs, float_epsilon, float_relative_epsilon); + } + private: float m_[9]; }; diff --git a/include/fire_engine/math/mat4.hpp b/include/fire_engine/math/mat4.hpp index 4bcae79c..4c0c3763 100644 --- a/include/fire_engine/math/mat4.hpp +++ b/include/fire_engine/math/mat4.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -115,13 +116,18 @@ class Mat4 return *this == rhs; } + // Approximate equality, component by component, through the ONE scalar authority + // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit + // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart + // than 1e-9, exactly as it always did), and both arguments mean both terms. An invalid + // tolerance — negative, NaN or infinite — makes the comparison FALSE rather than being + // reinterpreted. NaNs compare unequal now, which is the defect this replaced. [[nodiscard]] - constexpr bool approxEqual(const Mat4& rhs, float eps = float_epsilon) const noexcept + constexpr bool approxEqual(const Mat4& rhs, float eps, float relativeEps) const noexcept { for (int i = 0; i < 16; ++i) { - const float diff = m_[i] - rhs.m_[i]; - if (diff > eps || diff < -eps) + if (!almostEqual(m_[i], rhs.m_[i], eps, relativeEps)) { return false; } @@ -129,6 +135,19 @@ class Mat4 return true; } + [[nodiscard]] + constexpr bool approxEqual(const Mat4& rhs, float eps) const noexcept + { + // ABSOLUTE ONLY — a stated tolerance is the whole answer. + return approxEqual(rhs, eps, 0.0f); + } + + [[nodiscard]] + constexpr bool approxEqual(const Mat4& rhs) const noexcept + { + return approxEqual(rhs, float_epsilon, float_relative_epsilon); + } + [[nodiscard]] static Mat4 rotateY(float rad) noexcept { diff --git a/include/fire_engine/math/quaternion.hpp b/include/fire_engine/math/quaternion.hpp index e423ced2..4497a628 100644 --- a/include/fire_engine/math/quaternion.hpp +++ b/include/fire_engine/math/quaternion.hpp @@ -4,6 +4,7 @@ #include #include +#include #include namespace fire_engine @@ -101,15 +102,32 @@ class Quaternion return *this == rhs; } + // Approximate equality, component by component, through the ONE scalar authority + // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit + // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart + // than 1e-9, exactly as it always did), and both arguments mean both terms. An invalid + // tolerance — negative, NaN or infinite — makes the comparison FALSE rather than being + // reinterpreted. NaNs compare unequal now, which is the defect this replaced. [[nodiscard]] - constexpr bool approxEqual(const Quaternion& rhs, float eps = float_epsilon) const noexcept - { - const float dx = x_ - rhs.x_; - const float dy = y_ - rhs.y_; - const float dz = z_ - rhs.z_; - const float dw = w_ - rhs.w_; - return dx <= eps && dx >= -eps && dy <= eps && dy >= -eps && dz <= eps && dz >= -eps && - dw <= eps && dw >= -eps; + constexpr bool approxEqual(const Quaternion& rhs, float eps, float relativeEps) const noexcept + { + return almostEqual(x_, rhs.x_, eps, relativeEps) && + almostEqual(y_, rhs.y_, eps, relativeEps) && + almostEqual(z_, rhs.z_, eps, relativeEps) && + almostEqual(w_, rhs.w_, eps, relativeEps); + } + + [[nodiscard]] + constexpr bool approxEqual(const Quaternion& rhs, float eps) const noexcept + { + // ABSOLUTE ONLY — a stated tolerance is the whole answer. + return approxEqual(rhs, eps, 0.0f); + } + + [[nodiscard]] + constexpr bool approxEqual(const Quaternion& rhs) const noexcept + { + return approxEqual(rhs, float_epsilon, float_relative_epsilon); } [[nodiscard]] diff --git a/include/fire_engine/math/scalar.hpp b/include/fire_engine/math/scalar.hpp new file mode 100644 index 00000000..0ff9fea0 --- /dev/null +++ b/include/fire_engine/math/scalar.hpp @@ -0,0 +1,118 @@ +#pragma once + +#include + +#include + +// The ONE approximate scalar comparison, and the only thing in this library that decides whether +// two numbers are close enough (tier-0 review, finding 2). +// +// It exists because the pattern it replaces was wrong in a way that hid exactly the failures tests +// are for. Every vector, matrix and quaternion `approxEqual` computed +// +// const float diff = a - b; +// if (diff > eps || diff < -eps) { return false; } +// +// and BOTH comparisons are false when `diff` is NaN, so a NaN compared approximately equal to +// everything, itself included. Two equal infinities subtract to NaN and passed the same way — which +// reads as correct until you notice that opposite infinities do too. A test asserting a transform +// stayed finite could pass on a transform that was entirely NaN. + +namespace fire_engine +{ + +// The RELATIVE term's default. Deliberately its own constant rather than `float_epsilon`: that one +// is the absolute tolerance, and a single number serving as absolute tolerance, relative tolerance +// and degeneracy threshold is three unrelated policies wearing one name. +// +// 1e-6 is roughly ten times float's 1.19e-7 epsilon — close enough to machine precision to reject +// real error, loose enough to absorb the last couple of bits after a few operations. +inline constexpr float float_relative_epsilon = 1.0e-6f; + +// THREE OVERLOADS, because the tolerance a caller SPELLS OUT must mean what it says. +// +// almostEqual(a, b) default absolute AND default relative +// almostEqual(a, b, absolute) absolute only — exactly the historical comparison +// almostEqual(a, b, absolute, relative) both, stated +// +// A single function with a defaulted relative term would have made `almostEqual(a, b, 1e-9f)` admit +// a 1e-7 difference at magnitude 1.0: an explicit tolerance silently overridden by an implicit one. +// Overload selection decides the policy instead, so ordinary comparisons get the large-magnitude +// fix while an explicit tolerance remains the whole answer. +// +// TOLERANCES ARE VALIDATED, not reinterpreted. A negative, NaN or infinite tolerance is a caller +// defect — a misconfigured constant, an uninitialised field, a division that went wrong — and the +// comparison refuses it by returning FALSE, even for operands that are equal. Treating a negative +// as "a term that cannot be satisfied" is mathematically defensible and diagnostically useless: it +// converts a broken configuration into a different valid policy and lets the run continue. The +// check happens BEFORE the `a == b` shortcut precisely so equal operands cannot hide it. +// +// The comparison itself, once the tolerances are known good: +// +// 1. `a == b`, so identical values are equal whatever they are — including two equal infinities, +// which is the case the old subtraction form got right only by accident (inf - inf is NaN, and +// NaN passed its test). This also makes +0.0 and -0.0 equal, which is correct: they are the +// same number. +// 2. Any remaining non-finite OPERAND is not equal. After step 1 that means a NaN anywhere, or two +// infinities that differ, or an infinity against any finite value. NaN is never equal to +// anything, itself included — the property the old form inverted. +// 3. Finite values compare against the larger of the absolute and relative tolerances. The +// absolute term is what makes values near zero comparable at all, where relative error is +// meaningless; the relative term is what keeps large values comparable without a bespoke +// tolerance per call site. +// +// DOUBLE intermediates, and not for accuracy: `a - b` in float overflows to infinity for values +// near FLT_MAX of opposite sign, so the difference itself could be non-finite for two perfectly +// finite inputs, and the comparison would be answering about a number neither caller mentioned. In +// double the subtraction of any two floats is exact. +[[nodiscard]] constexpr bool almostEqual(float a, float b, float absoluteTolerance, + float relativeTolerance) noexcept +{ + // `std::isfinite`, not `(x - x) == 0`. The subtraction form classifies correctly but computes + // `inf - inf` on the way, which is an INVALID operation: it raises FE_INVALID and would trap + // where floating-point exceptions are enabled. A comparison must not alter exception state to + // answer a question about its arguments. (C++23 made these constexpr, so the function stays + // usable in constant expressions.) + const auto validTolerance = [](float tolerance) + { return std::isfinite(tolerance) && tolerance >= 0.0f; }; + if (!validTolerance(absoluteTolerance) || !validTolerance(relativeTolerance)) + { + return false; + } + + if (a == b) + { + return true; + } + if (!std::isfinite(a) || !std::isfinite(b)) + { + return false; + } + + const double difference = static_cast(a) - static_cast(b); + const double magnitude = difference < 0.0 ? -difference : difference; + if (magnitude <= static_cast(absoluteTolerance)) + { + return true; + } + const double scaleA = a < 0.0f ? -static_cast(a) : static_cast(a); + const double scaleB = b < 0.0f ? -static_cast(b) : static_cast(b); + const double scale = scaleA > scaleB ? scaleA : scaleB; + return magnitude <= static_cast(relativeTolerance) * scale; +} + +// ABSOLUTE ONLY. `almostEqual(a, b, 1e-9f)` rejects anything more than 1e-9 apart, at any +// magnitude, which is what it looks like it does. +[[nodiscard]] constexpr bool almostEqual(float a, float b, float absoluteTolerance) noexcept +{ + return almostEqual(a, b, absoluteTolerance, 0.0f); +} + +// The ordinary comparison: both defaults, so values near zero and values near FLT_MAX are each +// compared by the term that means something at their scale. +[[nodiscard]] constexpr bool almostEqual(float a, float b) noexcept +{ + return almostEqual(a, b, float_epsilon, float_relative_epsilon); +} + +} // namespace fire_engine diff --git a/include/fire_engine/math/vec_base.hpp b/include/fire_engine/math/vec_base.hpp index 4f75629c..4507da2e 100644 --- a/include/fire_engine/math/vec_base.hpp +++ b/include/fire_engine/math/vec_base.hpp @@ -4,6 +4,7 @@ #include #include +#include namespace fire_engine { @@ -158,13 +159,18 @@ class VecBase return self() == rhs; } + // Approximate equality, component by component, through the ONE scalar authority + // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit + // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart + // than 1e-9, exactly as it always did), and both arguments mean both terms. An invalid + // tolerance — negative, NaN or infinite — makes the comparison FALSE rather than being + // reinterpreted. NaNs compare unequal now, which is the defect this replaced. [[nodiscard]] - constexpr bool approxEqual(const Derived& rhs, float eps = float_epsilon) const noexcept + constexpr bool approxEqual(const Derived& rhs, float eps, float relativeEps) const noexcept { for (std::size_t i = 0; i < N; ++i) { - const float diff = data_[i] - rhs.data_[i]; - if (diff > eps || diff < -eps) + if (!almostEqual(data_[i], rhs.data_[i], eps, relativeEps)) { return false; } @@ -172,6 +178,19 @@ class VecBase return true; } + [[nodiscard]] + constexpr bool approxEqual(const Derived& rhs, float eps) const noexcept + { + // ABSOLUTE ONLY — a stated tolerance is the whole answer. + return approxEqual(rhs, eps, 0.0f); + } + + [[nodiscard]] + constexpr bool approxEqual(const Derived& rhs) const noexcept + { + return approxEqual(rhs, float_epsilon, float_relative_epsilon); + } + protected: float data_[N]{}; diff --git a/tests/math/test_scalar.cpp b/tests/math/test_scalar.cpp new file mode 100644 index 00000000..e4e0d435 --- /dev/null +++ b/tests/math/test_scalar.cpp @@ -0,0 +1,150 @@ +#include + +#include + +#include +#include + +#include +#include +#include +#include + +using namespace fire_engine; + +namespace +{ + +constexpr float kNaN = std::numeric_limits::quiet_NaN(); +constexpr float kInf = std::numeric_limits::infinity(); +constexpr float kMax = std::numeric_limits::max(); + +} // namespace + +TEST_CASE("a NaN is not approximately equal to anything, itself included", "[Scalar]") +{ + // THE defect this authority exists for (tier-0 finding 2). The previous form asked + // `diff > eps || diff < -eps`, and both comparisons are false for a NaN difference, so every + // approximate comparison in the engine answered "equal" for a NaN — including a test asserting + // that a transform had stayed finite. + CHECK_FALSE(almostEqual(kNaN, kNaN)); + CHECK_FALSE(almostEqual(kNaN, 0.0f)); + CHECK_FALSE(almostEqual(0.0f, kNaN)); + CHECK_FALSE(almostEqual(kNaN, kInf)); + // And no tolerance, however wide, may rescue it: NaN is not a value that is nearly something. + CHECK_FALSE(almostEqual(kNaN, 1.0f, 1.0e30f, 1.0e30f)); +} + +TEST_CASE("equal infinities are equal and opposite ones are not", "[Scalar]") +{ + // Both fell out of the old subtraction identically (inf - inf is NaN, and NaN passed), so the + // right answer and the wrong one were indistinguishable. Here they are decided before any + // arithmetic happens: `a == b` is the first question. + CHECK(almostEqual(kInf, kInf)); + CHECK(almostEqual(-kInf, -kInf)); + CHECK_FALSE(almostEqual(kInf, -kInf)); + // An infinity is not NEARLY a finite value, however large that value is. + CHECK_FALSE(almostEqual(kInf, kMax)); + CHECK_FALSE(almostEqual(-kInf, -kMax)); +} + +TEST_CASE("the difference is computed where it cannot overflow", "[Scalar]") +{ + // `kMax - (-kMax)` is +inf in float. Computed that way the comparison would be reasoning about + // a number neither caller passed; in double the subtraction of two floats is exact, so these + // answer about the values themselves. + CHECK_FALSE(almostEqual(kMax, -kMax)); + CHECK_FALSE(almostEqual(-kMax, kMax)); + // Two values a single ULP apart near FLT_MAX: far beyond any absolute tolerance, but well + // inside the relative one — which is the case the relative term exists for. + const float nextDown = std::nextafter(kMax, 0.0f); + CHECK(almostEqual(kMax, nextDown)); + CHECK_FALSE(almostEqual(kMax, nextDown, float_epsilon, 0.0f)); +} + +TEST_CASE("the tolerances are absolute and relative, in that order", "[Scalar]") +{ + // The FIRST argument keeps its historical meaning. A caller that wrote `approxEqual(rhs, 0.1f)` + // meant "within 0.1" and still does; reinterpreting it as 10% would have silently loosened + // every existing call site by orders of magnitude at large values. + CHECK(almostEqual(1.0f, 1.05f, 0.1f, 0.0f)); + CHECK_FALSE(almostEqual(1.0f, 1.5f, 0.1f, 0.0f)); + + // Near zero the relative term is meaningless (any two small values are relatively far apart), + // which is exactly where the absolute term carries the comparison. + CHECK(almostEqual(0.0f, 1.0e-9f)); + CHECK_FALSE(almostEqual(0.0f, 1.0e-9f, 0.0f, float_relative_epsilon)); + + // And far from zero the absolute term is meaningless, which is where the relative one does. + CHECK(almostEqual(1.0e6f, 1.0e6f + 0.5f)); + CHECK_FALSE(almostEqual(1.0e6f, 1.0e6f + 0.5f, float_epsilon, 0.0f)); +} + +TEST_CASE("+0.0 and -0.0 are the same number", "[Scalar]") +{ + // They compare equal here because they ARE equal numerically; a caller wanting to tell the two + // bit patterns apart wants something this function has never claimed to be. + CHECK(almostEqual(0.0f, -0.0f, 0.0f, 0.0f)); +} + +TEST_CASE("every vector, matrix and quaternion comparison inherits the NaN rule", "[Scalar]") +{ + // The authority is only worth having if nothing bypasses it. One NaN component per type, each + // of which used to compare equal to itself. + const Vec3 nanVec{kNaN, 0.0f, 0.0f}; + CHECK_FALSE(nanVec.approxEqual(nanVec)); + CHECK_FALSE(nanVec.approxEqual(Vec3{0.0f, 0.0f, 0.0f})); + + Mat3 nanMat3 = Mat3::identity(); + nanMat3[0, 0] = kNaN; + CHECK_FALSE(nanMat3.approxEqual(nanMat3)); + + Mat4 nanMat4 = Mat4::identity(); + nanMat4[2, 3] = kNaN; + CHECK_FALSE(nanMat4.approxEqual(nanMat4)); + + const Quaternion nanQuat{kNaN, 0.0f, 0.0f, 1.0f}; + CHECK_FALSE(nanQuat.approxEqual(nanQuat)); + + // Finite values still behave, so the fix is not simply "everything is unequal now". + CHECK(Vec3{1.0f, 2.0f, 3.0f}.approxEqual(Vec3{1.0f, 2.0f, 3.0f})); + CHECK(Mat4::identity().approxEqual(Mat4::identity())); +} + +TEST_CASE("an explicit tolerance is the whole answer", "[Scalar]") +{ + // The overload set exists for this. A caller who writes a tolerance means it: the two-argument + // form carries both defaults, the three-argument form is ABSOLUTE ONLY, and neither can be + // quietly loosened by the other's default. + CHECK(almostEqual(1.0f, 1.0f + 1.0e-7f)); // default relative admits it + CHECK_FALSE(almostEqual(1.0f, 1.0f + 1.0e-7f, 1.0e-9f)); // an explicit 1e-9 does not + CHECK(almostEqual(1.0f, 1.0f + 1.0e-7f, 1.0e-9f, 1.0e-6f)); + + // At large magnitudes the default form still works where an absolute tolerance cannot. + const float large = 1.0e7f; + CHECK(almostEqual(large, large + 1.0f)); + CHECK_FALSE(almostEqual(large, large + 1.0f, float_epsilon)); +} + +TEST_CASE("an invalid tolerance is refused, not reinterpreted", "[Scalar]") +{ + // A negative, NaN or infinite tolerance is a caller defect — a bad constant, an uninitialised + // field, a division that went wrong. Answering "false" makes it visible at the first + // comparison; treating a negative as an unsatisfiable term would silently convert a broken + // configuration into a stricter policy and let the run continue looking healthy. + CHECK_FALSE(almostEqual(1.0f, 1.0f, -1.0f)); + CHECK_FALSE(almostEqual(1.0f, 1.0f, 1.0e-6f, -1.0f)); + CHECK_FALSE(almostEqual(1.0f, 1.0f, kNaN)); + CHECK_FALSE(almostEqual(1.0f, 1.0f, 1.0e-6f, kNaN)); + CHECK_FALSE(almostEqual(1.0f, 1.0f, kInf)); + CHECK_FALSE(almostEqual(1.0f, 1.0f, 1.0e-6f, kInf)); + + // EQUAL OPERANDS TOO, which is why the validation runs before the `a == b` shortcut: the one + // case most likely to be exercised by a smoke test is the one that would hide the defect. + CHECK_FALSE(almostEqual(2.5f, 2.5f, -0.0001f)); + CHECK_FALSE(almostEqual(kInf, kInf, -1.0f)); + + // And the aggregates inherit the refusal rather than validating separately. + CHECK_FALSE(Vec3{1.0f, 2.0f, 3.0f}.approxEqual(Vec3{1.0f, 2.0f, 3.0f}, -1.0f)); + CHECK_FALSE(Mat4::identity().approxEqual(Mat4::identity(), kNaN)); +} From 58e283c508698c463ffacc2cba2aa1dcfdd2b6ef Mon Sep 17 00:00:00 2001 From: Nick Newson Date: Sun, 13 Sep 2026 15:33:29 +0100 Subject: [PATCH 2/5] Inversion answers whether it worked, not what it wishes were true MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Mat3::inverse() failed by returning the zero matrix — both a legitimate value and an error report, so a caller could not tell "singular" from "the answer is zero" without asking again. Its threshold was ABSOLUTE (|det| <= 1e-12), which rejects transforms that are perfectly well conditioned merely for being small: a uniform scale of 1e-5 has determinant 1e-15 and an exact inverse of 1e5, and the old test called it singular. VDPM had already worked around that with its own scale-invariant predicate, then called the absolute-threshold inverse anyway — so `coneUsable` reported true while `cameraObj` was a zero matrix times a camera offset. Two invertibility decisions, disagreeing, in adjacent lines. tryInverse() returns std::optional and is the only one. It normalises the matrix by its largest absolute component before taking the determinant, so the question asked is conditioning rather than size; it computes in double; it refuses non-finite input rather than producing a matrix of NaNs that every later operation spreads silently. Two things the optional means that a determinant test alone does not. An invalid tolerance is refused before the matrix is examined — `magnitude > tolerance` is TRUE for a zero determinant against a negative threshold, so a singular matrix would otherwise be accepted and then divided by its own zero. And every element is checked finite and within float's range in double BEFORE any conversion: a uniform scale of 1e-39 is perfectly conditioned and its inverse is 1e39, which float cannot hold, so converting would hand back an engaged optional full of infinities. The optional means "a representable inverse exists" or it means nothing. determinant() now returns double, and VDPM is why. A reflected instance with a tiny uniform scale has a usable inverse and a determinant that underflows float to -0.0f — which is not less than zero, so the facing sign says "winding preserved" for a transform that reverses it, and the cone culls the side that should be visible. The threshold for the underflow is lower than it looks: 1e-5 cubed is fine, 1e-16 cubed is not. VDPM's conditioning POLICY is reproduced exactly rather than inherited. tryInverse normalises by largest component while the old predicate used σ_max, and the two differ by (σ_max/maxComponent)³ ∈ [1, 27]; passing a constant would have relaxed the policy by up to 27× (or 27000× with the permissive default), and a blanket conservative constant would have REJECTED shapes the old test accepted. Scaling the tolerance by that exact ratio gives |det| > 1e-6·σ_max³ again — one decision, same answer, both sides of the boundary pinned by tests. The three physics inversions — the D block and Schur complement inside SpatialMatrix::inverse(), and invertDof's 3-DOF path — have a genuine invariant, so they go through a physics-local invertInvariant(): assert, and in release a logged reason and abort. An assert alone would leave *inverse dereferencing an empty optional in exactly the situation where the solver state is already corrupt. The helper stays in physics deliberately; in math/ it would recreate the unconditional inversion API this removes. Note for anyone reading the goldens: Determinism.GoldenHash builds a floor and three rigid boxes and never constructs an articulation, so it does not exercise these paths. The articulation suite is the relevant coverage. --- include/fire_engine/math/mat3.hpp | 151 +++++++++++++++++++---- include/fire_engine/physics/spatial.hpp | 45 ++++++- src/graphics/vdpm.cpp | 52 ++++++-- src/physics/articulation.cpp | 4 +- tests/graphics/test_vdpm.cpp | 57 +++++++++ tests/math/test_mat3.cpp | 153 ++++++++++++++++++++++-- 6 files changed, 415 insertions(+), 47 deletions(-) diff --git a/include/fire_engine/math/mat3.hpp b/include/fire_engine/math/mat3.hpp index cec2b39c..92217105 100644 --- a/include/fire_engine/math/mat3.hpp +++ b/include/fire_engine/math/mat3.hpp @@ -1,5 +1,9 @@ #pragma once +#include +#include +#include + #include #include #include @@ -7,6 +11,15 @@ namespace fire_engine { +// The conditioning threshold `tryInverse` applies to the NORMALISED determinant — its own named +// constant, because it is a policy about invertibility and not a comparison tolerance. +// +// 1e-9 sits about two orders BELOW float's epsilon (1.19e-7), so it is a permissive default: it +// rejects only what is numerically hopeless and leaves callers with a stricter requirement to say +// so. `makeVdpmViewParams` is one — its cone predicate wants a shape bound far tighter than this, +// and passes its own. +inline constexpr float kInverseConditionTolerance = 1.0e-9f; + // Column-major 3x3 matrix, mirroring Mat4's `[row, col]` accessor and storage // (`m_[col * 3 + row]`). Used for rotation matrices and (inverse) inertia tensors // in the rigid-body solver; kept minimal — only what the physics needs. @@ -102,42 +115,128 @@ class Mat3 return r; } + // DOUBLE, and the return type is the point rather than the intermediates. + // + // A float determinant of a tiny transform underflows: a uniform scale of 1e-5 has determinant + // 1e-15, and a tiny REFLECTED one lands at -0.0f — which reads as non-negative, so a caller + // deriving orientation from `det >= 0` concludes the transform preserves winding when it + // reverses it. VDPM does exactly that to fold a reflection into its cone facing, and the wrong + // answer there enables cone culling with an inverted facing sign: geometry culled from the side + // that should be visible. The magnitude is a conditioning question and the SIGN is an + // orientation question, and double keeps both answerable for transforms a float determinant + // cannot represent at all. [[nodiscard]] - constexpr float determinant() const noexcept + constexpr double determinant() const noexcept { - const float a = m_[0], b = m_[3], c = m_[6]; // row 0 - const float d = m_[1], e = m_[4], f = m_[7]; // row 1 - const float g = m_[2], h = m_[5], i = m_[8]; // row 2 + const double a = m_[0], b = m_[3], c = m_[6]; // row 0 + const double d = m_[1], e = m_[4], f = m_[7]; // row 1 + const double g = m_[2], h = m_[5], i = m_[8]; // row 2 return a * (e * i - f * h) + b * (f * g - d * i) + c * (d * h - e * g); } - // Inverse via the adjugate / determinant. Returns the zero matrix when (near-)singular - // (|det| <= eps) so callers can detect it and fall back rather than propagate NaNs. + // THE inverse, and the only one: there is no `inverse()` returning a zero matrix any more. + // + // The old signature failed by returning `Mat3{}`, which is both a legitimate value and an error + // report, so a caller could not tell "singular" from "the answer is zero" without checking + // again — and its threshold was ABSOLUTE (`|det| <= 1e-12`), which rejects transforms that are + // perfectly well conditioned merely for being small. A uniform scale of 1e-5 has determinant + // 1e-15 and an exact inverse of 1e5; the old test called it singular and handed back zeros. + // VDPM had already worked around this with its own scale-invariant predicate, then called the + // absolute-threshold inverse anyway and got the zero matrix while its own test said "usable". + // + // SCALE-INVARIANT, so the question asked is conditioning rather than size: the matrix is + // normalised by its largest absolute component before the determinant is taken, which puts a + // uniformly scaled transform and its unit-scale twin on exactly the same footing. `tolerance` + // is therefore a RELATIVE threshold on that normalised determinant, not a magnitude in the + // caller's units. + // + // Non-finite input is rejected rather than propagated: an inverse built from a NaN is a matrix + // of NaNs that every later operation quietly spreads. [[nodiscard]] - constexpr Mat3 inverse(float eps = 1.0e-12f) const noexcept + std::optional tryInverse(float tolerance = kInverseConditionTolerance) const noexcept { - const float a = m_[0], b = m_[3], c = m_[6]; // row 0 - const float d = m_[1], e = m_[4], f = m_[7]; // row 1 - const float g = m_[2], h = m_[5], i = m_[8]; // row 2 - const float A = e * i - f * h; - const float B = f * g - d * i; - const float C = d * h - e * g; - const float det = a * A + b * B + c * C; - if (det <= eps && det >= -eps) + // THE TOLERANCE IS VALIDATED FIRST, and not as a formality: `magnitude > tolerance` is TRUE + // for a zero determinant against a negative tolerance, so a singular matrix would be + // accepted and then divided by its own zero. A negative, NaN or infinite threshold is a + // caller defect, and the answer to it is "no inverse", never "every matrix is invertible". + if (!std::isfinite(tolerance) || tolerance < 0.0f) + { + return std::nullopt; + } + + double scale = 0.0; + for (const float value : m_) { - return Mat3{}; + if (!std::isfinite(value)) + { + return std::nullopt; + } + const double magnitude = + value < 0.0f ? -static_cast(value) : static_cast(value); + scale = magnitude > scale ? magnitude : scale; } - const float s = 1.0f / det; + if (scale == 0.0) + { + return std::nullopt; // the zero matrix: singular, and the one case scaling cannot help + } + + const double inverseScale = 1.0 / scale; + const double a = m_[0] * inverseScale, b = m_[3] * inverseScale, c = m_[6] * inverseScale; + const double d = m_[1] * inverseScale, e = m_[4] * inverseScale, f = m_[7] * inverseScale; + const double g = m_[2] * inverseScale, h = m_[5] * inverseScale, i = m_[8] * inverseScale; + const double A = e * i - f * h; + const double B = f * g - d * i; + const double C = d * h - e * g; + const double normalisedDet = a * A + b * B + c * C; + const double magnitude = normalisedDet < 0.0 ? -normalisedDet : normalisedDet; + if (!(magnitude > static_cast(tolerance))) + { + return std::nullopt; // `!(x > t)` so a NaN determinant is a refusal, not an acceptance + } + + // The normalised inverse, scaled back: inv(s·M) = inv(M)/s. + const double s = 1.0 / (normalisedDet * scale); + const double elements[9]{ + A * s, + B * s, + C * s, + (c * h - b * i) * s, + (a * i - c * g) * s, + (b * g - a * h) * s, + (b * f - c * e) * s, + (c * d - a * f) * s, + (a * e - b * d) * s, + }; + + // WELL CONDITIONED IS NOT THE SAME AS REPRESENTABLE. A uniform scale of 1e-39 is perfectly + // conditioned — its normalised determinant is 1 — and its inverse is 1e39, which float + // cannot hold. Converting anyway would hand back an engaged optional full of infinities + // (or, for values above float's range, an out-of-range conversion), so the caller would + // believe it had an inverse and propagate garbage. Every element is checked in double + // BEFORE any conversion happens, and the answer is "no usable inverse" instead. + // + // This is what makes the optional mean "a representable inverse exists", which is the only + // claim a caller can act on. + constexpr double kFloatMax = static_cast(std::numeric_limits::max()); + for (const double element : elements) + { + const double elementMagnitude = element < 0.0 ? -element : element; + if (!std::isfinite(element) || elementMagnitude > kFloatMax) + { + return std::nullopt; + } + } + Mat3 r; - r[0, 0] = A * s; - r[0, 1] = (c * h - b * i) * s; - r[0, 2] = (b * f - c * e) * s; - r[1, 0] = B * s; - r[1, 1] = (a * i - c * g) * s; - r[1, 2] = (c * d - a * f) * s; - r[2, 0] = C * s; - r[2, 1] = (b * g - a * h) * s; - r[2, 2] = (a * e - b * d) * s; + r[0, 0] = static_cast(elements[0]); + r[1, 0] = static_cast(elements[1]); + r[2, 0] = static_cast(elements[2]); + r[0, 1] = static_cast(elements[3]); + r[1, 1] = static_cast(elements[4]); + r[2, 1] = static_cast(elements[5]); + r[0, 2] = static_cast(elements[6]); + r[1, 2] = static_cast(elements[7]); + r[2, 2] = static_cast(elements[8]); return r; } diff --git a/include/fire_engine/physics/spatial.hpp b/include/fire_engine/physics/spatial.hpp index 6cf77921..efdb14ad 100644 --- a/include/fire_engine/physics/spatial.hpp +++ b/include/fire_engine/physics/spatial.hpp @@ -1,5 +1,10 @@ #pragma once +#include +#include +#include + +#include #include #include #include @@ -7,6 +12,41 @@ namespace fire_engine { +namespace physics_detail +{ + +// Invert a matrix the CALLER'S INVARIANT says must be invertible, and fail loudly if it is not. +// +// `Mat3::tryInverse` hands back an optional precisely so nobody can ignore a singular matrix. Two +// physics sites genuinely know better: an articulated-body inertia's `d` block is `m·1` plus +// positive-definite folded terms, and a 3-DOF joint's `D = SᵀU` is positive definite by +// construction. Neither can be singular unless the factorization upstream has already gone wrong. +// +// So the optional is checked, not assumed. An `assert` alone would leave `*inverse` dereferencing +// an empty optional in release — undefined behaviour in exactly the situation where the simulation +// state is already corrupt — so the release path is a defined stop: a logged reason, then +// `std::abort`. A wrong answer propagated through a solver is worse than a halt, because it comes +// back as a ragdoll that explodes three seconds later with nothing to point at. +// +// DELIBERATELY PHYSICS-LOCAL. Promoting this to `math/` would recreate the unconditional inversion +// API being removed: callers without an invariant would reach for it, and the optional would stop +// being the thing that makes them think. +[[nodiscard]] inline Mat3 invertInvariant(const Mat3& m, const char* what) +{ + const std::optional inverse = m.tryInverse(); + assert(inverse.has_value() && "a matrix the solver guarantees invertible was singular"); + if (!inverse.has_value()) + { + log::error(log::category::physics, + "{} was singular or non-finite; the articulated-body factorization is corrupt", + what); + std::abort(); + } + return *inverse; +} + +} // namespace physics_detail + // Minimal rigid (proper) transform — a unit quaternion rotation plus a translation, // mapping a point p ↦ rotation·p + translation. Unlike scene::Transform it carries no // scale or cached matrices: articulation forward kinematics composes thousands of these @@ -151,9 +191,10 @@ struct SpatialMatrix [[nodiscard]] SpatialMatrix inverse() const noexcept { - const Mat3 di = d.inverse(); + const Mat3 di = physics_detail::invertInvariant(d, "spatial inertia D block"); const Mat3 dic = di * c; - const Mat3 si = (a - b * dic).inverse(); + const Mat3 si = + physics_detail::invertInvariant(a - b * dic, "spatial inertia Schur complement"); const Mat3 bdi = b * di; const Mat3 negSiBdi = (si * bdi) * -1.0f; return SpatialMatrix{si, negSiBdi, (dic * si) * -1.0f, di - dic * negSiBdi}; diff --git a/src/graphics/vdpm.cpp b/src/graphics/vdpm.cpp index f8ac2119..975cf080 100644 --- a/src/graphics/vdpm.cpp +++ b/src/graphics/vdpm.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -455,20 +456,55 @@ VdpmViewParams makeVdpmViewParams(const Mat4& world, const Vec3& cameraPos, floa // The cone predicate runs in OBJECT space (sign of dot(normal, viewDir) is // transform-invariant), so the camera is inverse-transformed once. A near-singular world has no // reliable inverse — the cone is then unusable (never cull, treat every split as a potential - // silhouette). The conditioning test |det|/σ_max³ ∈ [0,1] is a pure SHAPE measure - // (scale-invariant), unlike an absolute |det| which would wrongly reject a tiny-but-uniform - // instance. A reflection's determinant sign folds into the cone facing. - const float det = linear.determinant(); - const float sigmaMax = worldLengthScale; - const bool coneUsable = std::abs(det) > 1e-6f * sigmaMax * sigmaMax * sigmaMax; + // silhouette). + // + // ONE DECISION, and it is the inverse itself. This used to ask a scale-invariant conditioning + // question here (|det| > 1e-6·σ_max³) and then call an inverse that applied its own ABSOLUTE + // threshold, so a tiny uniform instance passed the local test and got a zero matrix back from + // the inverse: `coneUsable` said yes while `cameraObj` was nonsense. `tryInverse` is now + // scale-invariant itself, so asking it IS the conditioning test, and there is no second opinion + // to disagree with. + // VDPM'S OWN THRESHOLD, and EXACTLY the one it always had — consolidating the decision must not + // change the policy, in either direction. + // + // The old predicate was |det| > 1e-6·σ_max³, a pure shape measure. `tryInverse` normalises by + // the largest absolute COMPONENT, so its threshold is on |det|/maxComponent³. Passing a + // constant would have shifted the policy by up to (σ_max/maxComponent)³ ∈ [1, 27]: 1e-6 admits + // matrices up to 27× closer to singular than before, and a blanket 27e-6 REJECTS ones the old + // test accepted (diag(1, 1, 1e-5) among them). Neither is "the same decision, consolidated". + // + // Scaling the tolerance by that exact ratio reproduces the original inequality: + // |det|/maxC³ > 1e-6·σ_max³/maxC³ ⟺ |det| > 1e-6·σ_max³ + // so one call now answers what two used to, with the same answer. + double maxComponent = 0.0; + for (int row = 0; row < 3; ++row) + { + for (int col = 0; col < 3; ++col) + { + maxComponent = std::max(maxComponent, std::abs(static_cast(linear[row, col]))); + } + } + // A zero matrix has no ratio to compute and no inverse either; `tryInverse` refuses it on its + // own, so the default tolerance is a placeholder rather than a policy in that case. + const double shapeRatio = + maxComponent > 0.0 ? static_cast(worldLengthScale) / maxComponent : 1.0; + const auto coneTolerance = static_cast(1.0e-6 * shapeRatio * shapeRatio * shapeRatio); + const std::optional linearInverse = linear.tryInverse(coneTolerance); + const bool coneUsable = linearInverse.has_value(); const Vec3 worldTranslation{world[0, 3], world[1, 3], world[2, 3]}; + // The sign comes from a DOUBLE determinant, and that is not a detail. A reflected instance with + // a tiny uniform scale has a perfectly usable inverse and a determinant that underflows float + // to -0.0f — which compares `>= 0.0f` as true, so the facing sign says "winding preserved" for + // a transform that reverses it, and the cone then culls the visible side. + const double det = linear.determinant(); + VdpmViewParams p; p.worldLinear = linear; p.worldTranslationMinusCamera = worldTranslation - cameraPos; - p.cameraObj = coneUsable ? linear.inverse() * (cameraPos - worldTranslation) : Vec3{}; + p.cameraObj = coneUsable ? *linearInverse * (cameraPos - worldTranslation) : Vec3{}; p.worldLengthScale = worldLengthScale; - p.facingSign = det >= 0.0f ? 1.0f : -1.0f; + p.facingSign = det >= 0.0 ? 1.0f : -1.0f; p.projScaleY = projScaleY; p.halfViewport = viewportHeight * 0.5f; p.silhouetteBoost = silhouetteBoost; diff --git a/src/physics/articulation.cpp b/src/physics/articulation.cpp index c4b4f501..18974c32 100644 --- a/src/physics/articulation.cpp +++ b/src/physics/articulation.cpp @@ -29,7 +29,9 @@ namespace } else if (nd == 3) { - r = d.inverse(); + // Positive definite by construction (D = SᵀU for a spherical joint), so a failure here is a + // corrupt factorization rather than a case to fall back from — see `invertInvariant`. + r = physics_detail::invertInvariant(d, "spherical joint D block"); } return r; } diff --git a/tests/graphics/test_vdpm.cpp b/tests/graphics/test_vdpm.cpp index 0afffffe..29e7cb12 100644 --- a/tests/graphics/test_vdpm.cpp +++ b/tests/graphics/test_vdpm.cpp @@ -690,6 +690,63 @@ TEST_CASE( CHECK(p.worldLengthScale < 10.01f); // and tight — the bound is exact for this symmetric case } +TEST_CASE("makeVdpmViewParams: a tiny reflected instance keeps its facing sign", "[vdpm]") +{ + // Two defects meet in this one transform, and both used to produce a plausible cone. + // + // FIRST, the inverse. A uniform scale of 1e-16 is perfectly conditioned, but its determinant is + // 1e-48. The old `Mat3::inverse()` applied an ABSOLUTE threshold (|det| <= 1e-12) and returned + // the zero matrix, while this function's own scale-invariant predicate said the cone WAS usable + // — so `coneUsable` was true and `cameraObj` was the zero vector times a camera offset. Two + // invertibility decisions, disagreeing. There is one now, and it is `tryInverse` itself. + // + // SECOND, the sign. A reflection reverses winding, which the cone folds into `facingSign`. In + // float that determinant underflows to -0.0f, which is not less than zero, so the sign reads + // +1 — the cone then culls the side that should be visible. The determinant is computed in + // double for exactly this reason. + Mat4 world = Mat4::identity(); + world[0, 0] = 1.0e-16f; + world[1, 1] = -1.0e-16f; // the reflection + world[2, 2] = 1.0e-16f; + + const VdpmViewParams p = + makeVdpmViewParams(world, Vec3{0, 0, 10}, 1.0f, 1000.0f, 0.0f, false, 1.0f, 1.0f, 1.0f); + + CHECK(p.coneUsable); // well conditioned, however small + CHECK(p.facingSign == -1.0f); + // And the object-space camera is a real inverse-transformed point rather than the zero vector + // the old path handed back: at this scale the camera is astronomically far away in object + // space, which is what a 1e-16 transform means. + CHECK(p.cameraObj.magnitude() > 1.0e10f); + + // The CONDITIONING POLICY is VDPM's own, and is EXACTLY the one the consolidation replaced: + // |det| > 1e-6·σ_max³. For diag(1, 1, z) the determinant is z and σ_max is 1, so the predicate + // reduces to z > 1e-6 — reproduced rather than approximated, because the cone's failure mode is + // culling geometry from the side that should be visible and its cost when disabled is only + // refinement work. + Mat4 conditioned = Mat4::identity(); + conditioned[2, 2] = 1.0e-5f; // accepted before the consolidation, and still accepted + const VdpmViewParams loose = makeVdpmViewParams(conditioned, Vec3{0, 0, 10}, 1.0f, 1000.0f, + 0.0f, false, 1.0f, 1.0f, 1.0f); + CHECK(loose.coneUsable); + // Rejected before, and still rejected — `tryInverse`'s own 1e-9 default would have taken this, + // which is the silent relaxation the explicit tolerance exists to prevent. + conditioned[2, 2] = 1.0e-7f; + const VdpmViewParams tight = makeVdpmViewParams(conditioned, Vec3{0, 0, 10}, 1.0f, 1000.0f, + 0.0f, false, 1.0f, 1.0f, 1.0f); + CHECK_FALSE(tight.coneUsable); + CHECK(tight.cameraObj == Vec3{}); // and an unusable cone reports no object-space camera at all + + // The unreflected twin differs ONLY in the sign, which is the point: conditioning and + // orientation are separate questions about the same matrix. + Mat4 unreflected = world; + unreflected[1, 1] = 1.0e-16f; + const VdpmViewParams q = makeVdpmViewParams(unreflected, Vec3{0, 0, 10}, 1.0f, 1000.0f, 0.0f, + false, 1.0f, 1.0f, 1.0f); + CHECK(q.coneUsable); + CHECK(q.facingSign == 1.0f); +} + TEST_CASE("refineForView: a singular world transform never back-face-culls", "[vdpm]") { // A degenerate (zero-scale) world has no reliable inverse, so the cone is unusable: diff --git a/tests/math/test_mat3.cpp b/tests/math/test_mat3.cpp index 6895ee6f..f6b1aca9 100644 --- a/tests/math/test_mat3.cpp +++ b/tests/math/test_mat3.cpp @@ -5,6 +5,9 @@ #include #include +#include +#include + using fire_engine::Mat3; using fire_engine::Quaternion; using fire_engine::Vec3; @@ -65,21 +68,151 @@ TEST_CASE("Mat3.InverseTimesMatrixIsIdentity", "[Mat3]") { // A non-symmetric invertible matrix (det != 0). const Mat3 a = Mat3::fromColumns({1.0f, 2.0f, 3.0f}, {0.0f, 1.0f, 4.0f}, {5.0f, 6.0f, 0.0f}); - const Mat3 ai = a.inverse(); - CHECK((a * ai).approxEqual(Mat3::identity(), 1e-4f)); - CHECK((ai * a).approxEqual(Mat3::identity(), 1e-4f)); + const std::optional ai = a.tryInverse(); + REQUIRE(ai.has_value()); + CHECK((a * *ai).approxEqual(Mat3::identity(), 1e-4f)); + CHECK((*ai * a).approxEqual(Mat3::identity(), 1e-4f)); // Identity inverts to itself; a diagonal inverts component-wise. - CHECK(Mat3::identity().inverse().approxEqual(Mat3::identity(), 1e-6f)); + REQUIRE(Mat3::identity().tryInverse().has_value()); + CHECK(Mat3::identity().tryInverse()->approxEqual(Mat3::identity(), 1e-6f)); + REQUIRE(Mat3::diagonal({2.0f, 4.0f, 0.5f}).tryInverse().has_value()); CHECK(Mat3::diagonal({2.0f, 4.0f, 0.5f}) - .inverse() - .approxEqual(Mat3::diagonal({0.5f, 0.25f, 2.0f}), 1e-6f)); + .tryInverse() + ->approxEqual(Mat3::diagonal({0.5f, 0.25f, 2.0f}), 1e-6f)); // A symmetric positive-definite matrix (the shape of the joint's effective-mass K). const Mat3 k = Mat3::fromColumns({4.0f, 1.0f, 0.5f}, {1.0f, 3.0f, 0.2f}, {0.5f, 0.2f, 2.0f}); - CHECK((k * k.inverse()).approxEqual(Mat3::identity(), 1e-4f)); + REQUIRE(k.tryInverse().has_value()); + CHECK((k * *k.tryInverse()).approxEqual(Mat3::identity(), 1e-4f)); +} + +TEST_CASE("Mat3.TryInverseAcceptsTinyWellConditionedTransforms", "[Mat3]") +{ + // THE tier-0 finding. A uniform scale of 1e-5 is perfectly conditioned — its inverse is a + // uniform 1e5 — but its determinant is 1e-15, which the old absolute threshold (|det| <= 1e-12) + // called singular and answered with a zero matrix. Size is not conditioning. + const Mat3 tiny = Mat3::diagonal({1.0e-5f, 1.0e-5f, 1.0e-5f}); + const std::optional inverse = tiny.tryInverse(); + REQUIRE(inverse.has_value()); + CHECK(inverse->approxEqual(Mat3::diagonal({1.0e5f, 1.0e5f, 1.0e5f}), 1.0f)); + CHECK((tiny * *inverse).approxEqual(Mat3::identity(), 1e-4f)); + + // The same shape across several decades: the answer must not depend on absolute size at all. + for (const float scale : {1.0e-8f, 1.0e-4f, 1.0f, 1.0e4f, 1.0e8f}) + { + const Mat3 scaled = Mat3::diagonal({scale, scale * 2.0f, scale * 0.5f}); + const std::optional scaledInverse = scaled.tryInverse(); + REQUIRE(scaledInverse.has_value()); + CHECK((scaled * *scaledInverse).approxEqual(Mat3::identity(), 1e-3f)); + } +} + +TEST_CASE("Mat3.TryInverseRejectsIllConditionedMatricesAtEveryScale", "[Mat3]") +{ + // Conditioning is a SHAPE question, so a matrix that is nearly rank-deficient must be refused + // whatever its magnitude — the mirror of the case above, and the reason the threshold applies + // to a normalised determinant rather than a raw one. + for (const float scale : {1.0e-6f, 1.0f, 1.0e6f}) + { + // Column 2 is (column 0 + column 1) to within a part in 1e11, so the NORMALISED determinant + // is ~1e-11 — an order below the 1e-9 conditioning threshold, whatever `scale` is. + const Mat3 nearlySingular = Mat3::fromColumns({scale, 0.0f, 0.0f}, {0.0f, scale, 0.0f}, + {scale, scale, scale * 1.0e-11f}); + CHECK_FALSE(nearlySingular.tryInverse().has_value()); + + // Calibration, so the threshold is a documented number rather than a mystery: the same + // shape two decades better conditioned is ACCEPTED at every scale. + const Mat3 conditioned = Mat3::fromColumns({scale, 0.0f, 0.0f}, {0.0f, scale, 0.0f}, + {scale, scale, scale * 1.0e-7f}); + CHECK(conditioned.tryInverse().has_value()); + } + + // Exactly singular: column 2 = column 0 + column 1. + const Mat3 singular = + Mat3::fromColumns({1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {1.0f, 1.0f, 0.0f}); + CHECK_FALSE(singular.tryInverse().has_value()); - // Singular matrix (rank-deficient: column 2 = column 0 + column 1) → zero, not NaN. - const Mat3 s = Mat3::fromColumns({1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {1.0f, 1.0f, 0.0f}); - CHECK(s.inverse() == Mat3{}); + // The zero matrix is the one singular case scaling cannot normalise, and it is refused rather + // than dividing by its own zero scale. + CHECK_FALSE(Mat3{}.tryInverse().has_value()); +} + +TEST_CASE("Mat3.TryInverseAcceptsReflectionsAndKeepsTheirSign", "[Mat3]") +{ + // A reflection is INVERTIBLE — its determinant is negative, not small — so a magnitude test is + // what the conditioning question needs, and the sign belongs to orientation instead. + const Mat3 reflection = Mat3::diagonal({1.0f, -1.0f, 1.0f}); + const std::optional inverse = reflection.tryInverse(); + REQUIRE(inverse.has_value()); + CHECK((reflection * *inverse).approxEqual(Mat3::identity(), 1e-6f)); + CHECK(reflection.determinant() < 0.0); + + // And a tiny reflection keeps that sign where a FLOAT determinant cannot. The underflow needs a + // scale small enough that the CUBE leaves float's range: 1e-5 cubed is -1e-15, which float + // represents perfectly well, but 1e-16 cubed is -1e-48 and float has nothing below about + // 1.4e-45. The sign is then lost to -0.0f, which is NOT less than zero and compares `>= 0` as + // true — so a caller reading the sign concludes the transform preserves winding when it + // reverses it. VDPM reads exactly this sign to fold a reflection into its cone facing, and the + // wrong answer culls the side that should be visible. + const Mat3 tinyReflection = Mat3::diagonal({1.0e-16f, -1.0e-16f, 1.0e-16f}); + REQUIRE(tinyReflection.tryInverse().has_value()); // conditioning is scale-invariant: usable + CHECK(tinyReflection.determinant() < 0.0); // in double the sign survives + const float asFloat = static_cast(tinyReflection.determinant()); + CHECK(asFloat == 0.0f); // in float the magnitude is gone... + CHECK_FALSE(asFloat < 0.0f); // ...and with it the orientation VDPM needs +} + +TEST_CASE("Mat3.TryInverseRefusesAnInvalidTolerance", "[Mat3]") +{ + // `magnitude > tolerance` is TRUE for a zero determinant against a negative tolerance, so + // without validation a singular matrix would be accepted and then divided by its own zero + // determinant. The tolerance is checked before the matrix is even examined. + const Mat3 singular = + Mat3::fromColumns({1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {1.0f, 1.0f, 0.0f}); + CHECK_FALSE(singular.tryInverse(-1.0f).has_value()); + CHECK_FALSE(singular.tryInverse(std::numeric_limits::quiet_NaN()).has_value()); + CHECK_FALSE(singular.tryInverse(std::numeric_limits::infinity()).has_value()); + + // A perfectly invertible matrix is refused too: a bad threshold is a caller defect, and the + // answer to one is never "every matrix is invertible". + CHECK_FALSE(Mat3::identity().tryInverse(-1.0f).has_value()); + CHECK_FALSE(Mat3::identity().tryInverse(std::numeric_limits::quiet_NaN()).has_value()); +} + +TEST_CASE("Mat3.TryInverseRefusesAnInverseFloatCannotHold", "[Mat3]") +{ + // WELL CONDITIONED IS NOT REPRESENTABLE. A uniform scale of 1e-39 has a normalised determinant + // of exactly 1 — it could not be better conditioned — and an inverse of 1e39, which is past + // float's 3.4e38. Converting anyway yields infinities inside an ENGAGED optional, so the caller + // believes it holds an inverse and propagates them. The optional means "a representable inverse + // exists" or it means nothing. + const Mat3 unrepresentable = Mat3::diagonal({1.0e-39f, 1.0e-39f, 1.0e-39f}); + CHECK_FALSE(unrepresentable.tryInverse().has_value()); + + // The neighbouring scale that IS representable still works, so the refusal is about the answer + // rather than about smallness. + const Mat3 representable = Mat3::diagonal({1.0e-38f, 1.0e-38f, 1.0e-38f}); + const std::optional inverse = representable.tryInverse(); + REQUIRE(inverse.has_value()); + CHECK((*inverse)[0, 0] > 0.0f); + CHECK(std::isfinite((*inverse)[0, 0])); + + // A mixed case: one axis fine, one beyond range. Any single unrepresentable element is enough, + // because the caller gets one matrix and cannot use half of it. + const Mat3 mixed = Mat3::diagonal({1.0f, 1.0e-39f, 1.0f}); + CHECK_FALSE(mixed.tryInverse().has_value()); +} + +TEST_CASE("Mat3.TryInverseRefusesNonFiniteMatrices", "[Mat3]") +{ + // An inverse built from a NaN is a matrix of NaNs that every later operation spreads silently. + // Refusing is what lets a caller notice at the point of failure. + Mat3 withNaN = Mat3::identity(); + withNaN[1, 1] = std::numeric_limits::quiet_NaN(); + CHECK_FALSE(withNaN.tryInverse().has_value()); + + Mat3 withInf = Mat3::identity(); + withInf[0, 2] = std::numeric_limits::infinity(); + CHECK_FALSE(withInf.tryInverse().has_value()); } From 478f29e1862bb3fe7115b19ae1c9b40d79729104 Mon Sep 17 00:00:00 2001 From: Nick Newson Date: Sun, 13 Sep 2026 16:07:57 +0100 Subject: [PATCH 3/5] Equality is exact, not bitwise, and says so now MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit bitwiseEqual() was `return self() == rhs` — ordinary floating-point comparison wearing a name that promised bit patterns. It was wrong in both directions it could be: -0.0f and +0.0f have different bits and compare EQUAL, while a NaN has bits identical to itself and compares UNEQUAL. It also duplicated operator== exactly, so it offered a choice between two spellings of one thing, one of which lied. Nothing outside math/ ever called it. Deleted rather than corrected: a method that means what the operator means earns its place only by being clearer, and this one was less clear. The comments claiming "strict bit-for-bit equality" are corrected on all four types — including Mat3, which made the claim without the method — and now state the two IEEE cases outright, plus what a caller would actually need if a determinism diagnostic ever wants real bit comparison (std::bit_cast, explicitly, at the point that needs it). The three tests were named BitwiseEqualMatchesOperator and asserted that the method agreed with the operator it forwarded to, which could not fail. They now assert the semantics the corrected comments claim: signed zeros compare equal, NaNs do not compare equal to themselves, and — for Quaternion — q and -q compare unequal despite being the same rotation, which is the seam the rotation redesign in phase 2 has to close. --- docs/onboarding.md | 14 ++++++++----- docs/review-order.md | 2 +- include/fire_engine/math/mat3.hpp | 6 +++++- include/fire_engine/math/mat4.hpp | 13 +++++------- include/fire_engine/math/quaternion.hpp | 18 ++++++++--------- include/fire_engine/math/vec_base.hpp | 13 +++++------- tests/math/test_mat4.cpp | 20 +++++++++++++++--- tests/math/test_quaternion.cpp | 17 +++++++++++++--- tests/math/test_vec3.cpp | 27 ++++++++++++++++++------- 9 files changed, 84 insertions(+), 46 deletions(-) diff --git a/docs/onboarding.md b/docs/onboarding.md index 4a90948e..d6cbc9fc 100644 --- a/docs/onboarding.md +++ b/docs/onboarding.md @@ -84,13 +84,17 @@ Start here because these classes are small, heavily tested, and used everywhere. - `Vec2`, `Vec3`, `Vec4`: numeric vector types with constexpr arithmetic and component accessors. `Vec3` also provides operations used by lighting, transforms, normals, and physics response. `magnitude()` / `normalise()` call `std::sqrt` and are intentionally - *not* `constexpr` (sqrt only became constexpr in C++26). `operator==` is strict bit - equality — use `approxEqual(rhs, eps)` for tolerance-based comparison (or `bitwiseEqual` - if you want to name the bit-identity intent explicitly). Vec3 ↔ Vec4 conversion is - `explicit` in both directions to prevent silent w-component loss/gain. + *not* `constexpr` (sqrt only became constexpr in C++26). `operator==` is **exact + component-wise IEEE equality — not bitwise**: `-0.0f` equals `+0.0f` though their bits + differ, and a NaN equals nothing though its bits are identical to itself. Use + `approxEqual(rhs, eps)` for tolerance-based comparison. (There was a `bitwiseEqual` that + simply called `operator==`; it claimed semantics it did not have and is gone. If you ever + need real bit identity, say so with `std::bit_cast` at the point that needs it.) + Vec3 ↔ Vec4 conversion is `explicit` in both directions to prevent silent w-component + loss/gain. - `Mat4`: column-major transform/projection matrix type. Look at translation, rotation, scale, perspective, and look-at helpers. Renderer, scene traversal, skinning, and physics - transforms all depend on this behaving predictably. Same `approxEqual` / `bitwiseEqual` + transforms all depend on this behaving predictably. Same `operator==` / `approxEqual` convention as the vector types. - `Quaternion`: runtime rotation representation for scene transforms. glTF rotations round trip better through quaternions than Euler angles. Animation uses SLERP for rotation diff --git a/docs/review-order.md b/docs/review-order.md index b9fcb4c3..c2e30a87 100644 --- a/docs/review-order.md +++ b/docs/review-order.md @@ -51,7 +51,7 @@ Read these first when a change touches build configuration, CI, or local tooling | `math/constants.hpp` | Just π/epsilon constants — orient quickly. | | `core/node_component_layout.hpp` + `.cpp` | Tiny but load-bearing: the rule deciding which of a glTF node's contents (Animator / Mesh / Light / Camera) owns the engine node and which move to identity-transform children. Exists because a `Node` holds ONE component while a glTF node may carry several, and the previous implicit rule — attach order — silently destroyed lights (`emplace` / `emplace` over an already-attached `Light`, no warning). Precedence is by what cannot move: an Animator must stay on the animated node, or its children stop following the animation. `materializeNodeComponentLayout` applies the plan to a real node — creating the identity children and returning the target for each payload — so the loader's attach sites consume a node rather than deciding placement from the current variant; that is what keeps the rule and the code from drifting, and it is Vulkan-free so CI verifies the actual topology. Exhaustively tested in `tests/core/test_node_component_layout.cpp`. | | `math/vec_base.hpp` | CRTP base for the vec types; compound-assign are primitives, binary ops delegate. | -| `math/vec2.hpp` / `vec3.hpp` / `vec4.hpp` | `Vec3` is the workhorse. **`magnitude()`/`normalise()` are deliberately NOT constexpr** (sqrt). `operator==` is **strict bit equality** — use `approxEqual` for tolerance. Vec3↔Vec4 conversions are `explicit` both ways. | +| `math/vec2.hpp` / `vec3.hpp` / `vec4.hpp` | `Vec3` is the workhorse. **`magnitude()`/`normalise()` are deliberately NOT constexpr** (sqrt). `operator==` is **exact component-wise IEEE equality, NOT bitwise** (`-0.0f` equals `+0.0f` with different bits; a NaN equals nothing with identical bits) — use `approxEqual` for tolerance, and `std::bit_cast` at the call site if bit identity is ever genuinely what is wanted. The `bitwiseEqual` that used to sit beside it merely called `operator==` and is removed. Vec3↔Vec4 conversions are `explicit` both ways. | | `math/quaternion.hpp` | SLERP, `fromVectors`, Hamilton `operator*`, and `integrate(ω, dt)` (exponential-map orientation integration for the rigid-body solver). Used for all scene rotation; glTF round-trips through this. | | `math/mat3.hpp` | Column-major 3×3 (mirrors `Mat4`'s `[row,col]`): `fromQuaternion`, `diagonal`, `transpose`, `Mat3·Mat3` / `Mat3·Vec3`. Holds the world inverse inertia `R·diag(invI)·Rᵀ` in the physics solver. | | `math/mat4.hpp` | **Column-major.** Translation/rotation/scale/perspective/look-at. Everything downstream trusts this — verify the multiplication and handedness conventions. | diff --git a/include/fire_engine/math/mat3.hpp b/include/fire_engine/math/mat3.hpp index 92217105..8d49799a 100644 --- a/include/fire_engine/math/mat3.hpp +++ b/include/fire_engine/math/mat3.hpp @@ -313,7 +313,11 @@ class Mat3 return r; } - // Strict bit-for-bit equality — use approxEqual for tolerance. + // EXACT component-wise IEEE equality — not bitwise, despite what this used to claim. Two + // differences matter and both are the float `==` operator's, not ours: `-0.0f` equals `+0.0f` + // though their bit patterns differ, and a NaN equals nothing at all though its bit pattern is + // identical to itself. Use `approxEqual` when you want tolerance; if a determinism diagnostic + // ever needs REAL bit comparison, it has to say so with `std::bit_cast`. [[nodiscard]] constexpr bool operator==(const Mat3& rhs) const noexcept { diff --git a/include/fire_engine/math/mat4.hpp b/include/fire_engine/math/mat4.hpp index 4c0c3763..ddce1482 100644 --- a/include/fire_engine/math/mat4.hpp +++ b/include/fire_engine/math/mat4.hpp @@ -95,8 +95,11 @@ class Mat4 return Vec3{r.x(), r.y(), r.z()}; } - // Strict bit-for-bit equality. Two matrices that differ by a single ULP - // compare not-equal — use approxEqual when you want tolerance. + // EXACT component-wise IEEE equality — not bitwise, despite what this used to claim. Two + // differences matter and both are the float `==` operator's, not ours: `-0.0f` equals `+0.0f` + // though their bit patterns differ, and a NaN equals nothing at all though its bit pattern is + // identical to itself. Use `approxEqual` when you want tolerance; if a determinism diagnostic + // ever needs REAL bit comparison, it has to say so with `std::bit_cast`. [[nodiscard]] constexpr bool operator==(const Mat4& rhs) const noexcept { @@ -110,12 +113,6 @@ class Mat4 return true; } - [[nodiscard]] - constexpr bool bitwiseEqual(const Mat4& rhs) const noexcept - { - return *this == rhs; - } - // Approximate equality, component by component, through the ONE scalar authority // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart diff --git a/include/fire_engine/math/quaternion.hpp b/include/fire_engine/math/quaternion.hpp index 4497a628..68083a71 100644 --- a/include/fire_engine/math/quaternion.hpp +++ b/include/fire_engine/math/quaternion.hpp @@ -86,22 +86,20 @@ class Quaternion return {-x_, -y_, -z_, -w_}; } - // Strict bit-for-bit equality. Two quaternions that differ by a single ULP - // compare not-equal — use approxEqual when you want tolerance. Note: also - // strict in the sign of the imaginary parts, so q and -q (which represent - // the same rotation) compare not-equal. + // EXACT component-wise IEEE equality — not bitwise, despite what this used to claim. Two + // differences matter and both are the float `==` operator's, not ours: `-0.0f` equals `+0.0f` + // though their bit patterns differ, and a NaN equals nothing at all though its bit pattern is + // identical to itself. Use `approxEqual` when you want tolerance; if a determinism diagnostic + // ever needs REAL bit comparison, it has to say so with `std::bit_cast`. + // + // Note this is a COMPONENT comparison, so `q` and `-q` are unequal here although they are the + // same rotation: the rotation-aware question belongs to a rotation type, not to this one. [[nodiscard]] constexpr bool operator==(const Quaternion& rhs) const noexcept { return x_ == rhs.x_ && y_ == rhs.y_ && z_ == rhs.z_ && w_ == rhs.w_; } - [[nodiscard]] - constexpr bool bitwiseEqual(const Quaternion& rhs) const noexcept - { - return *this == rhs; - } - // Approximate equality, component by component, through the ONE scalar authority // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart diff --git a/include/fire_engine/math/vec_base.hpp b/include/fire_engine/math/vec_base.hpp index 4507da2e..cce298ba 100644 --- a/include/fire_engine/math/vec_base.hpp +++ b/include/fire_engine/math/vec_base.hpp @@ -138,8 +138,11 @@ class VecBase return self(); } - // Strict bit-for-bit equality. Two values that differ by a single ULP - // compare not-equal — use approxEqual when you want tolerance. + // EXACT component-wise IEEE equality — not bitwise, despite what this used to claim. Two + // differences matter and both are the float `==` operator's, not ours: `-0.0f` equals `+0.0f` + // though their bit patterns differ, and a NaN equals nothing at all though its bit pattern is + // identical to itself. Use `approxEqual` when you want tolerance; if a determinism diagnostic + // ever needs REAL bit comparison, it has to say so with `std::bit_cast`. [[nodiscard]] friend constexpr bool operator==(const Derived& lhs, const Derived& rhs) noexcept { @@ -153,12 +156,6 @@ class VecBase return true; } - [[nodiscard]] - constexpr bool bitwiseEqual(const Derived& rhs) const noexcept - { - return self() == rhs; - } - // Approximate equality, component by component, through the ONE scalar authority // (`math/scalar.hpp`), and in its three forms — no argument means both defaults, an explicit // tolerance means ABSOLUTE ONLY (so `approxEqual(rhs, 1e-9f)` rejects anything further apart diff --git a/tests/math/test_mat4.cpp b/tests/math/test_mat4.cpp index fc77f307..8fc890b3 100644 --- a/tests/math/test_mat4.cpp +++ b/tests/math/test_mat4.cpp @@ -9,6 +9,8 @@ #include #include +#include + using fire_engine::Mat4; using fire_engine::Vec3; using fire_engine::Vec4; @@ -189,13 +191,25 @@ TEST_CASE("Mat4Equality.ZeroMatrices", "[Mat4Equality]") CHECK(a == b); } -TEST_CASE("Mat4Equality.BitwiseEqualMatchesOperator", "[Mat4Equality]") +TEST_CASE("Mat4Equality.EqualityIsExactNotBitwise", "[Mat4Equality]") { Mat4 a = Mat4::identity(); Mat4 b = Mat4::identity(); Mat4 c = Mat4::scale(Vec3{2.0f, 1.0f, 1.0f}); - CHECK(a.bitwiseEqual(b)); - CHECK_FALSE(a.bitwiseEqual(c)); + CHECK(a == b); + CHECK_FALSE(a == c); + // Exact component equality, not bitwise: -0.0f equals +0.0f though the bits differ, and a NaN + // equals nothing though its bits are identical to themselves. `bitwiseEqual()` claimed the + // latter semantics and delivered the former, so it is gone. + Mat4 signedZero = Mat4::identity(); + signedZero[0, 3] = -0.0f; + Mat4 plusZero = Mat4::identity(); + plusZero[0, 3] = 0.0f; + CHECK(signedZero == plusZero); + + Mat4 withNaN = Mat4::identity(); + withNaN[1, 1] = std::numeric_limits::quiet_NaN(); + CHECK_FALSE(withNaN == withNaN); } TEST_CASE("Mat4Equality.ApproxEqualWithinTolerance", "[Mat4Equality]") diff --git a/tests/math/test_quaternion.cpp b/tests/math/test_quaternion.cpp index 191619cb..9fc84685 100644 --- a/tests/math/test_quaternion.cpp +++ b/tests/math/test_quaternion.cpp @@ -6,6 +6,8 @@ #include #include +#include + #include #include #include @@ -96,13 +98,22 @@ TEST_CASE("QuaternionEquality.DifferentValuesNotEqual", "[QuaternionEquality]") CHECK_FALSE(a == b); } -TEST_CASE("QuaternionEquality.BitwiseEqualMatchesOperator", "[QuaternionEquality]") +TEST_CASE("QuaternionEquality.EqualityIsExactNotBitwise", "[QuaternionEquality]") { Quaternion a{0.1f, 0.2f, 0.3f, 0.4f}; Quaternion b{0.1f, 0.2f, 0.3f, 0.4f}; Quaternion c{0.1f, 0.2f, 0.3f, 0.5f}; - CHECK(a.bitwiseEqual(b)); - CHECK_FALSE(a.bitwiseEqual(c)); + CHECK(a == b); + CHECK_FALSE(a == c); + // Exact component equality, not bitwise — see Vec3Equality.EqualityIsExactNotBitwise. + CHECK(Quaternion{0.0f, 0.0f, 0.0f, 1.0f} == Quaternion{-0.0f, -0.0f, -0.0f, 1.0f}); + const float nan = std::numeric_limits::quiet_NaN(); + const Quaternion withNaN{nan, 0.0f, 0.0f, 1.0f}; + CHECK_FALSE(withNaN == withNaN); + + // And a COMPONENT comparison, so `q` and `-q` differ although they are the same rotation. + const Quaternion q{0.0f, 0.0f, 0.7071068f, 0.7071068f}; + CHECK_FALSE(q == Quaternion{-q.x(), -q.y(), -q.z(), -q.w()}); } TEST_CASE("QuaternionEquality.ApproxEqualWithinTolerance", "[QuaternionEquality]") diff --git a/tests/math/test_vec3.cpp b/tests/math/test_vec3.cpp index 526cc812..bfba4c24 100644 --- a/tests/math/test_vec3.cpp +++ b/tests/math/test_vec3.cpp @@ -139,13 +139,26 @@ TEST_CASE("Vec3Equality.NegativeZeroEqualsPositiveZero", "[Vec3Equality]") CHECK(a == b); } -TEST_CASE("Vec3Equality.BitwiseEqualMatchesOperator", "[Vec3Equality]") -{ - Vec3 a{1.0f, 2.0f, 3.0f}; - Vec3 b{1.0f, 2.0f, 3.0f}; - Vec3 c{1.0f, 2.0f, 3.5f}; - CHECK(a.bitwiseEqual(b)); - CHECK_FALSE(a.bitwiseEqual(c)); +TEST_CASE("Vec3Equality.EqualityIsExactNotBitwise", "[Vec3Equality]") +{ + // `operator==` is EXACT COMPONENT-WISE IEEE equality. It was described as "bit-for-bit" and + // duplicated by a `bitwiseEqual()` that simply called it, which was wrong in both directions — + // so both the name and the duplicate are gone, and the two cases that make the distinction real + // are pinned here instead. + const Vec3 a{1.0f, 2.0f, 3.0f}; + const Vec3 b{1.0f, 2.0f, 3.0f}; + const Vec3 c{1.0f, 2.0f, 3.5f}; + CHECK(a == b); + CHECK_FALSE(a == c); + + // DIFFERENT BITS, EQUAL VALUES: -0.0f and +0.0f have different sign bits and are the same + // number, so a genuinely bitwise comparison would answer the opposite of this. + CHECK(Vec3{0.0f, 0.0f, 0.0f} == Vec3{-0.0f, -0.0f, -0.0f}); + + // IDENTICAL BITS, UNEQUAL VALUES: a NaN never equals itself, whatever its payload. + const float nan = std::numeric_limits::quiet_NaN(); + const Vec3 withNaN{nan, 2.0f, 3.0f}; + CHECK_FALSE(withNaN == withNaN); } TEST_CASE("Vec3Equality.ApproxEqualWithinTolerance", "[Vec3Equality]") From 0dc596ac9b94bbb3815b4d1879f30c4656b9a976 Mon Sep 17 00:00:00 2001 From: Nick Newson Date: Sun, 13 Sep 2026 20:07:20 +0100 Subject: [PATCH 4/5] Norms that survive the ends of the float range, at no cost to the middle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit sqrt(dot(v, v)) fails at both ends. Components above ~1.8e19 square to infinity, so a vector whose length is perfectly representable reports infinity and then normalises to zero; components below ~1e-22 square to zero, so a small but ordinary vector reports a length of zero and loses its direction entirely. Neither case announces itself as a failure — both return a plausible number. Both failures DO announce themselves in the sum, which comes back infinite, NaN, zero, or subnormal. So the naive sum is computed first and trusted whenever it is finite and NORMAL, and the scaled fallback runs only where there was no correct answer before it. That ordering is the whole design: - Ordinary vectors take the arithmetic the engine always used, bit for bit. Neither physics golden moves. A change about extreme values has no business perturbing a box stack, and this one does not. - The robust path is confined to inputs that were already wrong. Normal, not merely positive, and that distinction is the difference between covering the underflow region and covering only its endpoint. A sum that has gone subnormal is finite and nonzero and has already lost most of its precision: (3e-23, 3e-23, 0) sums to 2.8e-45, carries about two significant bits, and answered 5.29e-23 against a true 4.24e-23 — 24.8% high, from a fast path that believed itself safe because the sum looked usable. A finite subnormal sum does reach the scaled fallback in normalise() as well; it simply cannot produce a normalised result, because such a magnitude is necessarily below float_normalise_cutoff and the configured degenerate value comes back instead. normalise() divides ONCE per component wherever the length is representable. The first version of this commit divided twice — by the largest component, then by a scaled norm — which is an extra ulp per component and looked harmless. It delayed a settling box stack from step 169 to 425 on macOS and from 167 to 1309 on Linux, against a 600-step budget, so the sleep demo failed on Linux only and looked for all the world like a marginal test. It was not marginal: friction directions that no longer cancel keep feeding the contact solver. A solver notices an ulp; that is what a solver is. The two-step form survives only for the case that needs it — a vector of finite components whose norm exceeds float's range, where dividing by the (correctly infinite) length would answer zero and throw away an ordinary direction. Three answers from normalise, and the middle one is the one that matters: a non-finite component yields NaNs, VISIBLY invalid; a magnitude below float_normalise_cutoff yields the zero vector or the identity rotation, as it always has; anything else is normalised. Laundering a corrupt input into the degenerate answer would turn a NaN orientation into a confident "no rotation" and surface it three seconds later somewhere unrelated. float_normalise_cutoff is its own constant at the same numeric value float_epsilon had. One number was serving as absolute comparison tolerance and degeneracy threshold at once; the rename moves no behaviour and stops the next person changing one policy while meaning the other. Measured, 2M calls: Quaternion::normalise 2.39 -> 2.49 ms (parity), Vec3::magnitude 3.27 -> 5.59, Vec3::normalise 2.68 -> 4.07. The residual on Vec3 is the guard branch defeating vectorisation in a synthetic back-to-back loop; swapping isfinite for a bare comparison measured identically (5.67 vs 5.59), so the clearer form stayed. Engine-level physics workloads are unchanged. --- CLAUDE.md | 1 + docs/onboarding.md | 29 ++++ docs/review-order.md | 2 +- include/fire_engine/math/constants.hpp | 10 ++ include/fire_engine/math/quaternion.hpp | 115 ++++++++++++++- include/fire_engine/math/vec_base.hpp | 177 ++++++++++++++++++++++-- tests/math/test_quaternion.cpp | 53 +++++++ tests/math/test_vec3.cpp | 102 ++++++++++++++ 8 files changed, 470 insertions(+), 19 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 22059187..3c32aa9a 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -150,6 +150,7 @@ assets/ glTF samples + HDR skyboxes - **One declaration of every shared GPU data-layout limit** — the sizes and indices that a C++ block and a shader block must agree on (caster/light/joint/morph/emitter/kernel counts, the map-validity bits). Purely GLSL-side algorithm constants are NOT in scope and this mechanism does not own them: a compute workgroup size, a scan radix, a tap count with no C++ counterpart stays where it is used. `shaders/gpu_limits.glsl` is written in the subset that is valid GLSL *and* valid C++; shaders `#include` it and `graphics/gpu_limits.hpp` includes it inside a `shader_limits` namespace, re-exporting each value under its `k`-name. Add a shader-visible limit **there**, never as a literal on either side, and keep the file inside the common subset (no `constexpr`, `inline`, `namespace`, `static_cast`, unsigned suffixes — each breaks the *other* language, in files that never mention this one). The `gpu_limits_guard` CTest case sweeps `shaders/` for a re-declaration, requires each consumer to use the name rather than a literal, and requires each `k`-constant to be defined *as* the shared declaration, so C++ cannot drift back to hard-coded values behind green shader checks. - **A shadow family's recording and its uploaded validity are one value** — `ShadowMapValidity` (`graphics/shadow_map_validity.hpp`) is applied twice per frame in `Renderer::prepareShadowPlan`, in a fixed order, both from the COMPLETED view set: as ELIGIBILITY, deciding which families may be PREPARED at all (preparation resolves casters and stages hysteresis, so a family that will be neither recorded nor sampled must not be resolved); then as CONFIRMATION (`shadowMapValidityFromPlan`), derived from the finished plan and judged against the counts eligibility expected, which is what `uploadFrameLighting` writes to `LightUBO::shadowMapValidMask` for every sampling path in `shader.frag`. Never skip a family's recording without routing the decision through it — a skipped family's depth image holds an earlier frame's content, and sampling it produces no error, no crash, and shadows from a frame that is gone. - **The shadow pass decides in preparation and records from the plan** — `prepareShadowFrame` (`graphics/shadow_pass_prepare.hpp`) filters, resolves each caster's LOD per view, claims the diagnostic row and builds a `ShadowFramePlan`; `Shadows::recordPass` consumes that plan and nothing else (no draw spans, no view set, no resolver). Anything the pass rasterises with must live in the prepared view or draw: a value read at record time that the comparison never saw is a cached shadow map kept when it should have been re-rendered. +- **One comparison authority, and norms with a scaled fallback** — `math/scalar.hpp`'s `almostEqual` is the only thing that decides whether two floats are close, and every `approxEqual` on `Vec*`/`Mat3`/`Mat4`/`Quaternion` delegates to it. `a == b` first (equal infinities pass), any remaining non-finite operand unequal (a NaN is never equal to anything, itself included), then absolute and relative tolerances in `double`. THREE OVERLOADS so an explicit tolerance is never loosened by an implicit one — none means both defaults, one is ABSOLUTE ONLY, two are both stated — and an invalid tolerance (negative, NaN, infinite) returns false rather than being reinterpreted, checked before the `a == b` shortcut so equal operands cannot hide a bad constant. `magnitude()`/`normalise()` compute `sqrt(dot(v, v))` FIRST and fall back to a scaled form only when that sum is not finite and NORMAL: the fast path keeps ordinary results bit-identical (so the physics goldens do not move for a change about extreme values), and the fallback covers overflow, total underflow AND the subnormal region, where a finite positive sum has already lost its precision. `normalise` divides ONCE per component wherever the length is representable — dividing twice costs an ulp each, which is not nothing: it tripled a box stack's settle time and only one platform's test noticed. Its three answers are load-bearing: non-finite → NaNs (VISIBLY invalid, never laundered into a plausible value), below `float_normalise_cutoff` (its own constant, not `float_epsilon`) → zero vector or identity rotation, otherwise normalised. `Mat3::tryInverse()` is the only inversion: scale-invariant, `double`, and its `std::optional` means "a REPRESENTABLE inverse exists" — an engaged value is never a matrix of infinities. Physics sites with a genuine invariant use the physics-local fail-fast helper, never an unchecked dereference. - **A reused shadow map is a claim about the GPU, so it is only ever made after the submit** — `ShadowResidencyStore` (`graphics/shadow_pass_plan.hpp`) records what each physical view's depth image HOLDS, and `prepareShadowFrame` compares this frame's prepared content against it to mark each view `Reused` or `Recorded`. It is owned by `Shadows`, beside the images it describes: that is the whole invalidation story, and why there is no `invalidate()` to forget to call — recreating the images means reconstructing the object that owns both. Two rules the type enforces rather than its callers: only a `Recorded` view commits (a `Reused` one never touched its image, so committing its prepared work would replace the record of what the image holds with a description of a frame that wrote nothing), and an `Invalid` slot is left alone (nothing recorded means nothing overwrote the image, so its record is still true). The commit sits beside `shadowLodResolver_.commitFrame()` BETWEEN `submitFrame` and `presentFrame`, for the same reason: content adopted by a frame that was abandoned would claim an image holds pixels the GPU never drew, and committing after PRESENTATION would be worse still — raii `presentKHR` throws on an out-of-date swapchain, so a resize would skip the commit for a frame whose depth was already being rasterised — and it is `noexcept`, adopting by MOVE out of the plan (`ShadowFramePlan::takeRecorded`, with `static_assert`s pinning the no-throw moves), because on the far side of a submit there is no useful answer to a failed allocation. `RenderTunables::shadowResidencyReuseEnabled` (overlay: "Reuse unchanged shadow views") forces every engaged view to record; it is SCHEDULING, so it is an argument to the law and never part of the content descriptor — a frame recorded with reuse off commits as usual and is reusable the moment it is switched back on. Each SH-01 row carries the disposition it ended up with, because zero raster passes alone cannot separate "reused" from "never engaged". - **GPU data-layout discipline** — every CPU struct shared with a shader (UBO/SSBO) lives in `render/ubo.hpp` with `alignas` + `static_assert`s pinning its std140/std430 offsets and size. Preserve this: when you change a shader-visible struct, update both sides and keep the static_asserts — they are the only thing catching a silent host↔GPU layout mismatch. Mapped host-visible writes go through `graphics/mapped_buffer.hpp` `writeMapped` (a bounds-checked `std::span`), never a raw `void*`. **And a block bound by more than one shader is declared ONCE, in a shared `shaders/*.glsl` include** (`light_ubo.glsl` for `LightUBO`, `material.glsl` for the bindless `Materials` SSBO + `MaterialData`, `shadow_push.glsl` for the `ShadowPushConstants` push block), never hand-copied per shader: field offsets depend on every field before them, so a copy missing an inserted field misreads everything after it, with no validation error and no crash. That is how the sky came to be multiplied by a shadow matrix — `selfShadowViewProj` was added to the struct and `shader.frag`, not to `skybox.frag`, and the wrong value read 1.0 until a scene had two skinned self-shadow casters. The `shader_block_guards` CTest case (`cmake/check_shader_blocks.cmake`) fails on a re-declared block *and* on a shared include that stops declaring it. diff --git a/docs/onboarding.md b/docs/onboarding.md index d6cbc9fc..a85aa6b4 100644 --- a/docs/onboarding.md +++ b/docs/onboarding.md @@ -961,6 +961,35 @@ the same change — most have a test or guard that will catch you, but not all. one constructed value so they cannot drift), and the per-frame-ring buffer handles are carried for recording but EXCLUDED from the comparison, since identical content alternates handles every frame. +- **Every approximate comparison goes through `almostEqual`, and every norm has a scaled fallback.** + `math/scalar.hpp` is the one place that decides whether two floats are close: `a == b` first (so + equal infinities pass), non-finite operands unequal (so a NaN is never equal to anything, itself + included), then absolute and relative tolerances in `double`. The types' `approxEqual` delegate; + none of them re-implements the test. Three overloads, because an explicit tolerance must not be + loosened by an implicit one — no argument means both defaults, ONE argument is absolute only, two + are both stated — and an invalid tolerance (negative, NaN, infinite) makes the comparison false + rather than being reinterpreted as a policy. + + `magnitude()` and `normalise()` compute `sqrt(dot(v, v))` FIRST and fall back to a scaled form + only when that sum comes back zero, infinite or NaN — which is exactly when the naive computation + had no answer (components above ~1.8e19 square to infinity; below ~1e-22 they flush to zero). Two + consequences to preserve if you touch this. Ordinary vectors take the arithmetic the engine always + used, BIT FOR BIT, so the physics goldens do not move for a change about extreme values. And the + fallback normalises through ONE division wherever the length is representable: dividing twice + (by the largest component, then by a scaled norm) costs an ulp per component, which delayed a + settling box stack from step 169 to 425 here and 1309 on Linux against a 600-step budget. + + Three answers from `normalise`, and the middle one is the one people get wrong: a non-finite input + yields NaNs (**visibly invalid**), a magnitude below `float_normalise_cutoff` yields the zero + vector or the identity rotation (degenerate, as documented), and anything else is normalised — + including a vector whose LENGTH is unrepresentable but whose direction is ordinary. Laundering an + invalid input into the degenerate answer is what makes a corrupt orientation surface three seconds + later somewhere unrelated. + + If you do change the arithmetic here, check `ReplayIsBitIdentical` and `FreeFallMatchesClosedForm` + before re-baselining any golden: those separate "last-bit arithmetic changed" from "the physics + changed", and a settle-time probe separates both from "the solver now takes three times as long to + come to rest". - **What an image HOLDS is committed only after the submit, and only for a view that recorded.** `ShadowResidencyStore` (`graphics/shadow_pass_plan.hpp`) is the other operand of the disposition law: preparation compares this frame's prepared content against it, and a view whose content diff --git a/docs/review-order.md b/docs/review-order.md index c2e30a87..7e02591c 100644 --- a/docs/review-order.md +++ b/docs/review-order.md @@ -51,7 +51,7 @@ Read these first when a change touches build configuration, CI, or local tooling | `math/constants.hpp` | Just π/epsilon constants — orient quickly. | | `core/node_component_layout.hpp` + `.cpp` | Tiny but load-bearing: the rule deciding which of a glTF node's contents (Animator / Mesh / Light / Camera) owns the engine node and which move to identity-transform children. Exists because a `Node` holds ONE component while a glTF node may carry several, and the previous implicit rule — attach order — silently destroyed lights (`emplace` / `emplace` over an already-attached `Light`, no warning). Precedence is by what cannot move: an Animator must stay on the animated node, or its children stop following the animation. `materializeNodeComponentLayout` applies the plan to a real node — creating the identity children and returning the target for each payload — so the loader's attach sites consume a node rather than deciding placement from the current variant; that is what keeps the rule and the code from drifting, and it is Vulkan-free so CI verifies the actual topology. Exhaustively tested in `tests/core/test_node_component_layout.cpp`. | | `math/vec_base.hpp` | CRTP base for the vec types; compound-assign are primitives, binary ops delegate. | -| `math/vec2.hpp` / `vec3.hpp` / `vec4.hpp` | `Vec3` is the workhorse. **`magnitude()`/`normalise()` are deliberately NOT constexpr** (sqrt). `operator==` is **exact component-wise IEEE equality, NOT bitwise** (`-0.0f` equals `+0.0f` with different bits; a NaN equals nothing with identical bits) — use `approxEqual` for tolerance, and `std::bit_cast` at the call site if bit identity is ever genuinely what is wanted. The `bitwiseEqual` that used to sit beside it merely called `operator==` and is removed. Vec3↔Vec4 conversions are `explicit` both ways. | +| `math/vec2.hpp` / `vec3.hpp` / `vec4.hpp` | `Vec3` is the workhorse. **`magnitude()`/`normalise()` are deliberately NOT constexpr** (sqrt) and compute `sqrt(dot(v, v))` FIRST, falling back to a scaled form only when that sum is not finite and NORMAL. The fast path is what keeps ordinary results bit-identical to the pre-tier-0 arithmetic (neither physics golden moved); the fallback is what makes `(1e20, 1e20, 0)` stop reporting an infinite length and `(3e-23, 3e-23, 0)` stop answering 24.8% high — a sum that has gone SUBNORMAL is finite and positive and has already lost its precision, which is why "positive" is not the guard. `normalise` divides ONCE per component wherever the length is representable: dividing twice costs an ulp each and delayed a settling box stack from step 169 to 425 (macOS) and 167 to 1309 (Linux). Its three answers are non-finite → NaNs (visibly invalid), below **`float_normalise_cutoff`** (its own constant, NOT `float_epsilon`) → zero vector / identity rotation, otherwise normalised — including a vector whose length overflows but whose direction is ordinary. `operator==` is **exact component-wise IEEE equality, NOT bitwise** (`-0.0f` equals `+0.0f` with different bits; a NaN equals nothing with identical bits) — use `approxEqual` for tolerance, and `std::bit_cast` at the call site if bit identity is ever genuinely what is wanted. The `bitwiseEqual` that used to sit beside it merely called `operator==` and is removed. Vec3↔Vec4 conversions are `explicit` both ways. | | `math/quaternion.hpp` | SLERP, `fromVectors`, Hamilton `operator*`, and `integrate(ω, dt)` (exponential-map orientation integration for the rigid-body solver). Used for all scene rotation; glTF round-trips through this. | | `math/mat3.hpp` | Column-major 3×3 (mirrors `Mat4`'s `[row,col]`): `fromQuaternion`, `diagonal`, `transpose`, `Mat3·Mat3` / `Mat3·Vec3`. Holds the world inverse inertia `R·diag(invI)·Rᵀ` in the physics solver. | | `math/mat4.hpp` | **Column-major.** Translation/rotation/scale/perspective/look-at. Everything downstream trusts this — verify the multiplication and handedness conventions. | diff --git a/include/fire_engine/math/constants.hpp b/include/fire_engine/math/constants.hpp index 203a9dbb..67eaac33 100644 --- a/include/fire_engine/math/constants.hpp +++ b/include/fire_engine/math/constants.hpp @@ -8,6 +8,16 @@ inline constexpr float deg_to_rad = pi / 180.0f; inline constexpr float rad_to_deg = 180.0f / pi; inline constexpr float float_epsilon = 1e-8f; +// Below this magnitude a vector or quaternion has no reliable DIRECTION, so `normalise` returns its +// documented degenerate answer (the zero vector, or the identity rotation) rather than dividing. +// +// Its own constant, deliberately. `float_epsilon` was doing this job as well as being the absolute +// comparison tolerance, which meant one number carried two unrelated policies and neither could be +// tuned without disturbing the other. The VALUE is unchanged, so this commit moves no behaviour +// with the rename; the point is that the next person changing a comparison tolerance does not +// silently change what counts as a degenerate direction. +inline constexpr float float_normalise_cutoff = 1e-8f; + // Soft pitch clamp for first-person cameras. Just under π/2 (≈85.94°) — keeps // the lookAt basis well-conditioned at the poles by avoiding the degenerate // straight-up / straight-down case while still letting the player look near diff --git a/include/fire_engine/math/quaternion.hpp b/include/fire_engine/math/quaternion.hpp index 68083a71..5389741d 100644 --- a/include/fire_engine/math/quaternion.hpp +++ b/include/fire_engine/math/quaternion.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include #include @@ -147,20 +148,49 @@ class Quaternion } [[nodiscard]] + // Fast path first — see `VecBase::magnitude`. While the sum of squares is finite and positive + // this is the arithmetic the engine always did, bit for bit; the scaled form below runs only + // when that sum came back zero, infinite or NaN, which is precisely when it had no answer. float magnitude() const noexcept { - return std::sqrt(magnitudeSquared()); + const float sumOfSquares = magnitudeSquared(); + // `isfinite` rather than a bare `< infinity` comparison: the two classify identically here + // (a NaN fails every comparison, an infinity fails the bound, zero fails the first test) + // and they measured identically too, so the one that says what it means wins. + // NORMAL, not merely positive. A sum that has gone SUBNORMAL has already lost most of its + // precision without reaching zero: (3e-23, 3e-23, 0) sums to 2.8e-45, which carries about + // two significant bits, and `sqrt` of it answers 5.29e-23 against a true 4.24e-23 — a 24.8% + // error from a fast path that thought it was fine because the sum was finite and positive. + // Requiring the sum to be at least `float`'s smallest NORMAL value routes that whole region + // to the scaled form, where the components are rescaled before they are squared and no + // precision is lost at all. `isfinite` then rules out the top end; a NaN fails both. + if (sumOfSquares >= std::numeric_limits::min() && std::isfinite(sumOfSquares)) + { + return std::sqrt(sumOfSquares); + } + return scaledMagnitude(); } + // A rotation's norm is the number every unit-quaternion assumption rests on — rotate(), + // slerp(), toMat4() — so the same three answers as the vector types. A non-finite component + // yields a NaN quaternion (visibly invalid, never laundered into the identity), a magnitude + // below `float_normalise_cutoff` yields the identity rotation as it always has, and anything + // else is normalised. [[nodiscard]] static Quaternion normalise(const Quaternion& q) noexcept { - float len = q.magnitude(); - if (len < float_epsilon) + const float sumOfSquares = q.magnitudeSquared(); + // A NORMAL sum, for the subnormal-precision reason given on `magnitude`. + if (sumOfSquares >= std::numeric_limits::min() && std::isfinite(sumOfSquares)) { - return Quaternion::identity(); + const float length = std::sqrt(sumOfSquares); + if (length < float_normalise_cutoff) + { + return Quaternion::identity(); + } + return {q.x_ / length, q.y_ / length, q.z_ / length, q.w_ / length}; } - return {q.x_ / len, q.y_ / len, q.z_ / len, q.w_ / len}; + return scaledNormalise(q); } Quaternion& normalise() noexcept @@ -404,6 +434,81 @@ class Quaternion } private: + // THE ROBUST PATH, reached only when the sum of squares was zero, infinite or NaN. + [[nodiscard]] float scaledMagnitude() const noexcept + { + const float components[4]{x_, y_, z_, w_}; + float largest = 0.0f; + bool anyInfinite = false; + for (const float component : components) + { + if (std::isnan(component)) + { + return component; + } + if (std::isinf(component)) + { + anyInfinite = true; + continue; + } + const float componentMagnitude = std::fabs(component); + largest = componentMagnitude > largest ? componentMagnitude : largest; + } + if (anyInfinite) + { + return std::numeric_limits::infinity(); + } + if (largest == 0.0f) + { + return 0.0f; + } + float sumOfScaledSquares = 0.0f; + for (const float component : components) + { + const float scaled = component / largest; + sumOfScaledSquares += scaled * scaled; + } + return largest * std::sqrt(sumOfScaledSquares); + } + + [[nodiscard]] static Quaternion scaledNormalise(const Quaternion& q) noexcept + { + const float components[4]{q.x_, q.y_, q.z_, q.w_}; + float largest = 0.0f; + for (const float component : components) + { + if (!std::isfinite(component)) + { + const float nan = std::numeric_limits::quiet_NaN(); + return {nan, nan, nan, nan}; + } + const float componentMagnitude = std::fabs(component); + largest = componentMagnitude > largest ? componentMagnitude : largest; + } + if (largest == 0.0f) + { + return Quaternion::identity(); + } + float sumOfScaledSquares = 0.0f; + for (const float component : components) + { + const float scaled = component / largest; + sumOfScaledSquares += scaled * scaled; + } + const float scaledNorm = std::sqrt(sumOfScaledSquares); // in [1, 2] + const float length = largest * scaledNorm; + if (length < float_normalise_cutoff) + { + return Quaternion::identity(); + } + if (std::isfinite(length)) + { + return {q.x_ / length, q.y_ / length, q.z_ / length, q.w_ / length}; + } + return {(q.x_ / largest) / scaledNorm, (q.y_ / largest) / scaledNorm, + (q.z_ / largest) / scaledNorm, (q.w_ / largest) / scaledNorm}; + } + float x_{0.0f}; float y_{0.0f}; float z_{0.0f}; diff --git a/include/fire_engine/math/vec_base.hpp b/include/fire_engine/math/vec_base.hpp index cce298ba..d8d33892 100644 --- a/include/fire_engine/math/vec_base.hpp +++ b/include/fire_engine/math/vec_base.hpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -98,13 +99,51 @@ class VecBase return Derived::dotProduct(self(), rhs); } - // std::sqrt is not constexpr before C++26, so magnitude (and the - // normalise helpers that depend on it) cannot be constexpr either. + // A FAST PATH THAT IS THE OLD ARITHMETIC, and a scaled fallback for the cases where the old + // arithmetic was wrong. + // + // `sqrt(dot(v, v))` fails at both ends of float's range: components above ~1.8e19 square to + // infinity, so a vector whose norm is perfectly representable reports infinity; components + // below ~1e-22 square to zero, so a small but ordinary vector reports zero and then normalises + // away to nothing. Both failures ANNOUNCE THEMSELVES in the sum — it comes back infinite, NaN, + // or exactly zero — so the naive sum can be computed first and trusted whenever it is finite + // and positive, which is every vector a frame of this engine actually contains. + // + // That matters for two reasons beyond speed (the scaled form measured ~3x the cost of this + // one). It keeps ordinary results BIT-IDENTICAL to what the engine computed before, so the + // physics goldens do not move for a change that was supposed to be about extreme values. And + // it confines the robust path to inputs where there was no correct answer before it. + // + // Special values are decided rather than inherited: any NaN component makes the norm NaN + // (checked first, so a vector holding both a NaN and an infinity is NaN rather than depending + // on iteration order); any infinity with no NaN makes it infinite; an all-zero vector is zero. + // + // std::sqrt is not constexpr before C++26, so this and the normalise helpers cannot be either. [[nodiscard]] float magnitude() const noexcept { - return std::sqrt(magnitudeSquared()); + const float sumOfSquares = magnitudeSquared(); + // `isfinite` rather than a bare `< infinity` comparison: the two classify identically here + // (a NaN fails every comparison, an infinity fails the bound, zero fails the first test) + // and they measured identically too, so the one that says what it means wins. + // NORMAL, not merely positive. A sum that has gone SUBNORMAL has already lost most of its + // precision without reaching zero: (3e-23, 3e-23, 0) sums to 2.8e-45, which carries about + // two significant bits, and `sqrt` of it answers 5.29e-23 against a true 4.24e-23 — a 24.8% + // error from a fast path that thought it was fine because the sum was finite and positive. + // Requiring the sum to be at least `float`'s smallest NORMAL value routes that whole region + // to the scaled form, where the components are rescaled before they are squared and no + // precision is lost at all. `isfinite` then rules out the top end; a NaN fails both. + if (sumOfSquares >= std::numeric_limits::min() && std::isfinite(sumOfSquares)) + { + return std::sqrt(sumOfSquares); + } + return scaledMagnitude(); } + // The raw sum of squares, which OVERFLOWS where `magnitude()` does not — for large components + // it is infinity and for tiny ones zero. That is honest for what it is (a squared quantity has + // half the exponent range available to it), and it stays because comparisons of squared lengths + // are a legitimate and cheaper thing to want. Reach for `magnitude()` when the answer is the + // length itself. [[nodiscard]] constexpr float magnitudeSquared() const noexcept { float sum = 0.0f; @@ -115,21 +154,40 @@ class VecBase return sum; } + // Normalised through the same fast path, for the same reasons: while the sum of squares is + // finite and positive, this is exactly the division the engine did before — one rounding per + // component, bit-identical results. Rounding twice instead (dividing by the largest component + // and then by a scaled norm) costs an ulp per component, which sounds like nothing and delayed + // a settling box stack from step 169 to 425 on macOS and 1309 on Linux against a 600-step + // budget. A solver notices an ulp; that is what solvers are. + // + // Three answers, and each is deliberate: + // * a non-finite component yields a NaN vector — VISIBLY invalid. Returning the zero vector + // or some identity would launder a corrupt input into a plausible value. + // * a magnitude below `float_normalise_cutoff` yields the zero vector, as it always has: + // there is no direction to report. + // * anything else is normalised — including a vector whose LENGTH is unrepresentable but + // whose direction is ordinary, which is the case that needs the scaled form. [[nodiscard]] static Derived normalise(const Derived& v) noexcept { - float len = v.magnitude(); - if (len < float_epsilon) - { - return Derived{}; - } - - Derived result{v}; - for (std::size_t i = 0; i < N; ++i) + const float sumOfSquares = v.magnitudeSquared(); + // A NORMAL sum, for the subnormal-precision reason given on `magnitude`. + if (sumOfSquares >= std::numeric_limits::min() && std::isfinite(sumOfSquares)) { - result.data_[i] /= len; + const float length = std::sqrt(sumOfSquares); + if (length < float_normalise_cutoff) + { + return Derived{}; + } + Derived result{v}; + for (std::size_t i = 0; i < N; ++i) + { + result.data_[i] /= length; + } + return result; } - return result; + return scaledNormalise(v); } Derived& normalise() noexcept @@ -188,6 +246,99 @@ class VecBase return approxEqual(rhs, float_epsilon, float_relative_epsilon); } +protected: + // THE ROBUST PATH, reached only when the sum of squares was zero, infinite or NaN — i.e. when + // the naive computation had no answer to give. Kept out of line from the hot path above. + [[nodiscard]] float scaledMagnitude() const noexcept + { + float largest = 0.0f; + bool anyInfinite = false; + for (std::size_t i = 0; i < N; ++i) + { + if (std::isnan(data_[i])) + { + return data_[i]; // a NaN outranks an infinity elsewhere in the vector + } + if (std::isinf(data_[i])) + { + anyInfinite = true; + continue; + } + const float componentMagnitude = std::fabs(data_[i]); + largest = componentMagnitude > largest ? componentMagnitude : largest; + } + if (anyInfinite) + { + return std::numeric_limits::infinity(); + } + if (largest == 0.0f) + { + return 0.0f; + } + + float sumOfScaledSquares = 0.0f; + for (std::size_t i = 0; i < N; ++i) + { + const float scaled = data_[i] / largest; + sumOfScaledSquares += scaled * scaled; + } + return largest * std::sqrt(sumOfScaledSquares); + } + + [[nodiscard]] static Derived scaledNormalise(const Derived& v) noexcept + { + float largest = 0.0f; + for (std::size_t i = 0; i < N; ++i) + { + if (!std::isfinite(v.data_[i])) + { + Derived invalid{}; + for (std::size_t j = 0; j < N; ++j) + { + invalid.data_[j] = std::numeric_limits::quiet_NaN(); + } + return invalid; + } + const float componentMagnitude = std::fabs(v.data_[i]); + largest = componentMagnitude > largest ? componentMagnitude : largest; + } + if (largest == 0.0f) + { + return Derived{}; + } + + float sumOfScaledSquares = 0.0f; + for (std::size_t i = 0; i < N; ++i) + { + const float scaled = v.data_[i] / largest; + sumOfScaledSquares += scaled * scaled; + } + const float scaledNorm = std::sqrt(sumOfScaledSquares); // in [1, sqrt(N)] + const float length = largest * scaledNorm; + if (length < float_normalise_cutoff) + { + return Derived{}; + } + if (std::isfinite(length)) + { + Derived result{v}; + for (std::size_t i = 0; i < N; ++i) + { + result.data_[i] = v.data_[i] / length; + } + return result; + } + // The length itself is unrepresentable — a vector of finite components whose norm exceeds + // float's range. The direction is still well defined, and dividing by an infinite length + // would answer zero, so the two-step scaled form is the only way to keep it. + Derived result{v}; + for (std::size_t i = 0; i < N; ++i) + { + result.data_[i] = (v.data_[i] / largest) / scaledNorm; + } + return result; + } + protected: float data_[N]{}; diff --git a/tests/math/test_quaternion.cpp b/tests/math/test_quaternion.cpp index 9fc84685..0c170d54 100644 --- a/tests/math/test_quaternion.cpp +++ b/tests/math/test_quaternion.cpp @@ -459,3 +459,56 @@ TEST_CASE("Quaternion.IntegrateAdvancesOrientation", "[Quaternion]") CHECK(stepped.rotate(Vec3{1.0f, 0.0f, 0.0f}) .approxEqual(expected.rotate(Vec3{1.0f, 0.0f, 0.0f}), 1e-4f)); } + +TEST_CASE("Quaternion.NormIsRobustAtBothEndsOfTheRange", "[Quaternion]") +{ + // The same scaled norm as VecBase, and it matters more here: every unit-quaternion assumption + // in the engine — rotate(), slerp(), toMat4() — rests on this one number. + const Quaternion huge{1.0e20f, 1.0e20f, 0.0f, 0.0f}; + CHECK(huge.magnitudeSquared() == std::numeric_limits::infinity()); + CHECK(std::isfinite(huge.magnitude())); + CHECK(huge.magnitude() == Catch::Approx(1.41421356e20f).epsilon(1e-5)); + + const Quaternion tiny{1.0e-25f, 1.0e-25f, 0.0f, 0.0f}; + CHECK(tiny.magnitudeSquared() == 0.0f); + CHECK(tiny.magnitude() > 0.0f); + + // A quaternion whose length float cannot hold still has a direction, so normalisation works + // from the scaled components rather than dividing by an infinity. + const Quaternion vast{3.0e38f, 3.0e38f, 3.0e38f, 3.0e38f}; + CHECK(std::isinf(vast.magnitude())); + const Quaternion unit = Quaternion::normalise(vast); + CHECK(unit.magnitude() == Catch::Approx(1.0f).epsilon(1e-6)); + CHECK(unit.x() == Catch::Approx(0.5f).epsilon(1e-6)); +} + +TEST_CASE("Quaternion.NormIsAccurateWhereTheSumGoesSubnormal", "[Quaternion]") +{ + // The same hole as Vec3.MagnitudeIsAccurateWhereTheSumGoesSUBNORMAL: a sum that is positive and + // finite but SUBNORMAL has already lost most of its bits, so "finite and positive" is not a + // sufficient guard for the fast path. + const Quaternion tiny{3.0e-23f, 3.0e-23f, 0.0f, 0.0f}; + REQUIRE(tiny.magnitudeSquared() > 0.0f); + REQUIRE(tiny.magnitudeSquared() < std::numeric_limits::min()); + CHECK(tiny.magnitude() == Catch::Approx(4.2426407e-23f).epsilon(1e-6)); +} + +TEST_CASE("Quaternion.NormaliseDistinguishesDegenerateFromInvalid", "[Quaternion]") +{ + const float nan = std::numeric_limits::quiet_NaN(); + const float inf = std::numeric_limits::infinity(); + + // DEGENERATE keeps the documented answer: a quaternion with no length has no rotation, and the + // identity is the safe one to return. + CHECK(Quaternion::normalise(Quaternion{0.0f, 0.0f, 0.0f, 0.0f}) == Quaternion::identity()); + + // INVALID must stay visibly invalid. Returning the identity for a NaN input would turn a + // corrupt orientation into a confident "no rotation", which is exactly the kind of laundering + // that makes a physics bug surface three seconds later somewhere else. + const Quaternion fromNaN = Quaternion::normalise(Quaternion{nan, 0.0f, 0.0f, 1.0f}); + CHECK(std::isnan(fromNaN.w())); + CHECK_FALSE(fromNaN == Quaternion::identity()); + const Quaternion fromInf = Quaternion::normalise(Quaternion{0.0f, inf, 0.0f, 1.0f}); + CHECK(std::isnan(fromInf.w())); + CHECK_FALSE(fromInf == Quaternion::identity()); +} diff --git a/tests/math/test_vec3.cpp b/tests/math/test_vec3.cpp index bfba4c24..00d5c50b 100644 --- a/tests/math/test_vec3.cpp +++ b/tests/math/test_vec3.cpp @@ -822,3 +822,105 @@ TEST_CASE("Vec3Compound.PlusEqualsChained", "[Vec3Compound]") (a += {1.0f, 0.0f, 0.0f}) += {0.0f, 1.0f, 0.0f}; expectNear(a, 2.0f, 2.0f, 1.0f); } + +TEST_CASE("Vec3.MagnitudeSurvivesScalesThatOverflowTheSquares", "[Vec3]") +{ + // `sqrt(dot(v, v))` fails at both ends of float's range, and a robust norm is the fix for both. + // + // TOO LARGE: (1e20, 1e20, 0) squares to 1e40 per component, which overflows float — so the + // naive form reports an infinite length for a vector whose true length (1.41e20) is perfectly + // representable. Everything downstream then normalises to zero. + const Vec3 huge{1.0e20f, 1.0e20f, 0.0f}; + CHECK(huge.magnitudeSquared() == std::numeric_limits::infinity()); // the honest squared + CHECK(std::isfinite(huge.magnitude())); + CHECK(huge.magnitude() == Catch::Approx(1.41421356e20f).epsilon(1e-5)); + + // TOO SMALL: (1e-25, 1e-25, 0) squares to 1e-50, which flushes to zero — so the naive form + // reports a length of zero for a vector that is small but entirely ordinary, and the direction + // is lost rather than merely imprecise. + const Vec3 tiny{1.0e-25f, 1.0e-25f, 0.0f}; + CHECK(tiny.magnitudeSquared() == 0.0f); // again honest for a squared quantity + CHECK(tiny.magnitude() > 0.0f); + CHECK(tiny.magnitude() == Catch::Approx(1.41421356e-25f).epsilon(1e-5)); +} + +TEST_CASE("Vec3.MagnitudeIsAccurateWhereTheSumGoesSUBNORMAL", "[Vec3]") +{ + // The gap a "finite and positive" fast path leaves behind. (3e-23, 3e-23, 0) squares to 9e-46 + // per component and sums to 2.8e-45 — not zero, so it looks like a usable sum, but subnormal + // and carrying about two significant bits. `sqrt` of that answers 5.29e-23 against a true + // 4.24e-23: a 24.8% error from a path that believed itself safe. + // + // The fast path therefore requires a NORMAL sum, and everything below it is rescaled before + // squaring, which loses nothing. + const Vec3 subnormalSum{3.0e-23f, 3.0e-23f, 0.0f}; + REQUIRE(subnormalSum.magnitudeSquared() > 0.0f); // the sum survives... + REQUIRE(subnormalSum.magnitudeSquared() < + std::numeric_limits::min()); // ...but subnormal + CHECK(subnormalSum.magnitude() == Catch::Approx(4.2426407e-23f).epsilon(1e-6)); + + // Either side of the switch, so the boundary itself is covered rather than one point near it. + // 1e-19 squares to a normal sum and takes the fast path; the rest go subnormal and do not. + for (const float component : {1.0e-19f, 1.0e-20f, 1.0e-22f, 1.0e-23f, 1.0e-25f}) + { + const Vec3 v{component, component, 0.0f}; + const float expected = component * std::sqrt(2.0f); + CHECK(v.magnitude() == Catch::Approx(expected).epsilon(1e-5)); + } + + // NORMALISE IS A DIFFERENT QUESTION, and the answer here is the cutoff's, not the norm's. A + // finite subnormal sum DOES take the scaled fallback in `normalise` — the control flow is the + // same — but its magnitude is necessarily below `float_normalise_cutoff` (1e-8), since a + // subnormal sum of squares implies a magnitude of about 1e-19 at most. So the fallback runs and + // then returns the configured degenerate value anyway. The subnormal fix therefore changes what + // `magnitude()` answers and never what `normalise()` answers: the two share a fast path, not a + // policy. + CHECK(Vec3::normalise(Vec3{3.0e-23f, 3.0e-23f, 0.0f}) == Vec3{}); + CHECK(Vec3::normalise(Vec3{1.0e-19f, 1.0e-19f, 0.0f}) == Vec3{}); +} + +TEST_CASE("Vec3.MagnitudeDecidesItsSpecialValues", "[Vec3]") +{ + const float nan = std::numeric_limits::quiet_NaN(); + const float inf = std::numeric_limits::infinity(); + + CHECK(std::isnan(Vec3{nan, 1.0f, 2.0f}.magnitude())); + CHECK(std::isinf(Vec3{inf, 1.0f, 2.0f}.magnitude())); + CHECK(Vec3{-inf, 0.0f, 0.0f}.magnitude() > 0.0f); // magnitude is unsigned + // A NaN OUTRANKS an infinity, whichever order they appear in — otherwise the answer would + // depend on which component the loop reached first. + CHECK(std::isnan(Vec3{inf, nan, 0.0f}.magnitude())); + CHECK(std::isnan(Vec3{nan, inf, 0.0f}.magnitude())); + CHECK(Vec3{}.magnitude() == 0.0f); +} + +TEST_CASE("Vec3.NormaliseKeepsADirectionWhoseLengthIsUnrepresentable", "[Vec3]") +{ + // The reason normalisation works from the SCALED components rather than dividing by + // `magnitude()`: this vector's true length is about 5.2e38, past float's 3.4e38, so the length + // genuinely is infinity — while the direction is an ordinary diagonal. Dividing by that + // infinity yields the zero vector, silently turning a well-defined direction into nothing. + const Vec3 vast{3.0e38f, 3.0e38f, 3.0e38f}; + CHECK(std::isinf(vast.magnitude())); // correctly infinite: the LENGTH is not representable + const Vec3 direction = Vec3::normalise(vast); + const float expected = 1.0f / std::sqrt(3.0f); + CHECK(direction.approxEqual(Vec3{expected, expected, expected}, 1e-6f)); + CHECK(direction.magnitude() == Catch::Approx(1.0f).epsilon(1e-6)); +} + +TEST_CASE("Vec3.NormaliseReportsInvalidInputAsInvalid", "[Vec3]") +{ + // A non-finite input must NOT be laundered into a plausible value. Returning the zero vector + // here would look exactly like the documented degenerate answer for a zero-length vector, and + // the caller would carry a corrupt direction forward believing it was merely small. + const float nan = std::numeric_limits::quiet_NaN(); + const float inf = std::numeric_limits::infinity(); + const Vec3 fromNaN = Vec3::normalise(Vec3{nan, 1.0f, 0.0f}); + CHECK(std::isnan(fromNaN.x())); + const Vec3 fromInf = Vec3::normalise(Vec3{inf, 1.0f, 0.0f}); + CHECK(std::isnan(fromInf.x())); + + // While a genuinely degenerate input keeps the documented answer. + CHECK(Vec3::normalise(Vec3{}) == Vec3{}); + CHECK(Vec3::normalise(Vec3{1.0e-30f, 0.0f, 0.0f}) == Vec3{}); +} From dc7a2ede27258563740082644bf5a4d6254172f0 Mon Sep 17 00:00:00 2001 From: Nick Newson Date: Sun, 13 Sep 2026 20:10:22 +0100 Subject: [PATCH 5/5] Assert how long the stack takes to sleep, not merely that it does MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Demos.Sleep checked that everything was asleep by step 600. That deadline was generous enough never to flake and therefore generous enough to hide a regression: when the tier-0 norm work added an ulp of error per component, the settle time went from step 169 to 425 on macOS and the test still passed, with 175 steps to spare. Only Linux failed, because its trajectory was slower still — 1309 — and it looked exactly like a marginal test on a slow platform. The obvious response, widening the budget, would have buried the defect permanently. The step count is now the evidence. The test records the first post-impact step at which the whole island — the three-box stack AND the striker — is asleep, and bounds it at 260. That number is measured, not chosen for comfort: healthy runs settle at 167 (macOS/arm64) and 172 (Linux/x86_64), so 260 leaves better than 50% headroom over the slower platform while sitting decisively below the 425 a real regression produced. The gap between those two figures is what makes the assertion worth having; a bound inside it can distinguish them. It then steps a further 120 and re-checks, which separates "crossed the sleep threshold once" from "came to rest". A body oscillating around the threshold satisfies the first and fails the second. Widen the bound only with a measurement showing the engine legitimately got slower, and never to make a red test green — the comment in the test says so, because that is exactly the pressure this assertion will apply. The count is in ABSOLUTE simulation steps, and that is not incidental. Every figure here and in the branch history — 167, 172, 425, 1309 — is an absolute step index, so a loop counting from the end of the impact window would mean something quite different: 260 iterations starting after step 150 is a bound of 410, which is above the regression this test exists to catch. It reported settling at 19 where every measurement says 169. The loop therefore runs from kImpactWindowEnd + 1 to kSettleByStep and reports the absolute step. Also closes tier-0 phase 1 in the docs: docs/codereview.md records the resolution against findings 1, 2, 4 and 5 as five commits (this one included, since it exists only because of how the norms regression was caught), and docs/roadmap.md drops phase 1 entirely — it indexes open work, and this phase is no longer open. --- docs/codereview.md | 53 +++++++++++++++++++++++++ docs/onboarding.md | 7 +++- docs/review-order.md | 5 ++- docs/roadmap.md | 21 ++++++---- include/fire_engine/math/quaternion.hpp | 11 +++-- include/fire_engine/math/vec_base.hpp | 8 ++-- tests/physics/test_demos.cpp | 47 ++++++++++++++++++++-- 7 files changed, 130 insertions(+), 22 deletions(-) diff --git a/docs/codereview.md b/docs/codereview.md index 5f884af9..b896e43b 100644 --- a/docs/codereview.md +++ b/docs/codereview.md @@ -261,6 +261,59 @@ property/invariant tests for: The existing `Mat3` tests particularly need expansion: inversion is a critical operation, but the current suite does not cover scale invariance or the absolute-determinant failure above. +### Phase 1 resolution — ✅ landed (`math-correctness-foundation`, five commits) + +Findings 1, 2, 4 and 5 are cleared. Finding 3 (the rotation invariant) and findings 6–8 belong to +phases 2–3 and remain open. + +1. **One comparison authority** (finding 2). `math/scalar.hpp`'s `almostEqual` is the only thing that + decides whether two floats are close: `a == b` first (equal infinities pass), any remaining + non-finite operand unequal (a NaN is never equal to anything, itself included), then absolute and + relative tolerances in `double` so the difference itself cannot overflow. Every type's + `approxEqual` delegates. THREE OVERLOADS, because an explicit tolerance must not be loosened by + an implicit one — none means both defaults, one is absolute only, two are both stated. Invalid + tolerances (negative, NaN, infinite) return false rather than being reinterpreted, and the check + precedes the `a == b` shortcut so equal operands cannot hide a bad constant. +2. **`Mat3::tryInverse()`** (finding 1). Scale-invariant, `double` intermediates, `std::optional` + instead of a zero-matrix sentinel that was both a value and an error report. It refuses an + invalid tolerance before examining the matrix (`magnitude > tolerance` is TRUE for a zero + determinant against a negative threshold — a singular matrix accepted, then divided by its own + zero), and refuses an inverse `float` cannot represent (a `1e-39` uniform scale is perfectly + conditioned and its inverse is `1e39`, so an engaged optional would have held infinities). The + three physics invariant sites use a physics-local fail-fast helper — assert, then a logged reason + and `abort` in release, never an unchecked dereference. VDPM's two disagreeing invertibility + decisions became one, with its `|det| > 1e-6·σ_max³` policy reproduced exactly rather than + inherited; `determinant()` returns `double` so a tiny reflection's sign survives instead of + underflowing to `-0.0f` and inverting the cone facing. +3. **Equality is exact, not bitwise** (finding 5). `bitwiseEqual` was `return self() == rhs` — it + duplicated `operator==` and misdescribed it in both directions (`-0.0f` equals `+0.0f` with + different bits; a NaN equals nothing with identical bits). Deleted, comments corrected on all + four types, and the tests now assert those two IEEE cases plus the `q` vs `-q` seam that phase 2 + must close. +4. **Robust norms** (finding 4), and the shape of this one is the finding worth keeping. + `sqrt(dot(v, v))` is computed FIRST and trusted only when the sum is finite and **normal**; the + scaled form runs otherwise. That ordering keeps ordinary vectors bit-identical to the previous + arithmetic — **neither physics golden moved** — while covering both failure regions. Requiring a + normal sum rather than merely a positive one matters: `(3e-23, 3e-23, 0)` sums to `2.8e-45`, + which is finite, positive and carries about two significant bits, and answered 24.8% high. A + finite subnormal sum does reach the scaled fallback in `normalise` as well, but its magnitude is + necessarily below `float_normalise_cutoff`, so normalisation still returns the configured + degenerate value. `normalise` divides ONCE per component wherever the length is representable: + the first attempt divided twice, and that single extra ulp per component delayed a settling box + stack from step 169 to 425 on macOS and 167 to 1309 on Linux. +5. **A settle-time tripwire** (no finding — it exists because of how (4) was caught). Every endpoint + assertion in the suite passed on macOS while the broken version tripled the settle time; only + Linux failed, and only because its trajectory was slower still. `Demos.Sleep` now records the + first post-impact step at which the whole island sleeps and bounds it at 260 — measured from + healthy runs of 167 (macOS/arm64) and 172 (Linux/x86_64), decisively below the 425 a real + regression produced — then steps a further 120 to separate "crossed the threshold once" from + stable rest. + +The lesson worth carrying into phases 2 and 3: a green suite plus a dutifully re-baselined golden +looked exactly like success while the change was a regression. What separated them was measuring a +QUANTITY (settle steps) rather than asserting an endpoint, and treating a golden move as a question +rather than a chore. + ### Recommended implementation sequence #### Phase 1: correctness foundation diff --git a/docs/onboarding.md b/docs/onboarding.md index a85aa6b4..a962dadd 100644 --- a/docs/onboarding.md +++ b/docs/onboarding.md @@ -971,8 +971,11 @@ the same change — most have a test or guard that will catch you, but not all. rather than being reinterpreted as a policy. `magnitude()` and `normalise()` compute `sqrt(dot(v, v))` FIRST and fall back to a scaled form - only when that sum comes back zero, infinite or NaN — which is exactly when the naive computation - had no answer (components above ~1.8e19 square to infinity; below ~1e-22 they flush to zero). Two + whenever that sum is not finite and NORMAL — zero, SUBNORMAL, infinite or NaN — which is exactly + when the naive computation had no accurate answer (components above ~1.8e19 square to infinity; + below ~1e-22 they flush to zero; and in between, a subnormal sum is finite and positive while + carrying only a couple of significant bits, which is how `(3e-23, 3e-23, 0)` answered 24.8% high). + Two consequences to preserve if you touch this. Ordinary vectors take the arithmetic the engine always used, BIT FOR BIT, so the physics goldens do not move for a change about extreme values. And the fallback normalises through ONE division wherever the length is representable: dividing twice diff --git a/docs/review-order.md b/docs/review-order.md index 7e02591c..d7e57a8f 100644 --- a/docs/review-order.md +++ b/docs/review-order.md @@ -48,11 +48,12 @@ Read these first when a change touches build configuration, CI, or local tooling | File | Pay attention to | |---|---| -| `math/constants.hpp` | Just π/epsilon constants — orient quickly. | +| `math/constants.hpp` | π and the angle conversions, plus TWO tolerances that are deliberately separate policies: `float_epsilon` is the absolute comparison tolerance, and `float_normalise_cutoff` is the magnitude below which a vector or quaternion has no direction worth reporting. They hold the same value today and answer different questions; one number doing both jobs is what the tier-0 split fixed, so changing a comparison tolerance cannot silently redefine degeneracy. | +| `math/scalar.hpp` | **The one approximate comparison in the library** — every `approxEqual` on `Vec*`/`Mat3`/`Mat4`/`Quaternion` delegates here, and nothing re-implements the test. `a == b` first (so equal infinities pass), any remaining non-finite operand UNEQUAL (so a NaN is never equal to anything, itself included — the defect this replaced inverted exactly that), then absolute and relative tolerances compared in `double`, since `a - b` in float overflows near FLT_MAX and would answer about a number neither caller passed. THREE OVERLOADS, because an explicit tolerance must not be loosened by an implicit one: none means both defaults, one is ABSOLUTE ONLY, two are both stated. An invalid tolerance (negative, NaN, infinite) returns false rather than being reinterpreted as a policy, and that check runs BEFORE the `a == b` shortcut so equal operands cannot hide a bad constant. | | `core/node_component_layout.hpp` + `.cpp` | Tiny but load-bearing: the rule deciding which of a glTF node's contents (Animator / Mesh / Light / Camera) owns the engine node and which move to identity-transform children. Exists because a `Node` holds ONE component while a glTF node may carry several, and the previous implicit rule — attach order — silently destroyed lights (`emplace` / `emplace` over an already-attached `Light`, no warning). Precedence is by what cannot move: an Animator must stay on the animated node, or its children stop following the animation. `materializeNodeComponentLayout` applies the plan to a real node — creating the identity children and returning the target for each payload — so the loader's attach sites consume a node rather than deciding placement from the current variant; that is what keeps the rule and the code from drifting, and it is Vulkan-free so CI verifies the actual topology. Exhaustively tested in `tests/core/test_node_component_layout.cpp`. | | `math/vec_base.hpp` | CRTP base for the vec types; compound-assign are primitives, binary ops delegate. | | `math/vec2.hpp` / `vec3.hpp` / `vec4.hpp` | `Vec3` is the workhorse. **`magnitude()`/`normalise()` are deliberately NOT constexpr** (sqrt) and compute `sqrt(dot(v, v))` FIRST, falling back to a scaled form only when that sum is not finite and NORMAL. The fast path is what keeps ordinary results bit-identical to the pre-tier-0 arithmetic (neither physics golden moved); the fallback is what makes `(1e20, 1e20, 0)` stop reporting an infinite length and `(3e-23, 3e-23, 0)` stop answering 24.8% high — a sum that has gone SUBNORMAL is finite and positive and has already lost its precision, which is why "positive" is not the guard. `normalise` divides ONCE per component wherever the length is representable: dividing twice costs an ulp each and delayed a settling box stack from step 169 to 425 (macOS) and 167 to 1309 (Linux). Its three answers are non-finite → NaNs (visibly invalid), below **`float_normalise_cutoff`** (its own constant, NOT `float_epsilon`) → zero vector / identity rotation, otherwise normalised — including a vector whose length overflows but whose direction is ordinary. `operator==` is **exact component-wise IEEE equality, NOT bitwise** (`-0.0f` equals `+0.0f` with different bits; a NaN equals nothing with identical bits) — use `approxEqual` for tolerance, and `std::bit_cast` at the call site if bit identity is ever genuinely what is wanted. The `bitwiseEqual` that used to sit beside it merely called `operator==` and is removed. Vec3↔Vec4 conversions are `explicit` both ways. | -| `math/quaternion.hpp` | SLERP, `fromVectors`, Hamilton `operator*`, and `integrate(ω, dt)` (exponential-map orientation integration for the rigid-body solver). Used for all scene rotation; glTF round-trips through this. | +| `math/quaternion.hpp` | SLERP, `fromVectors`, Hamilton `operator*`, and `integrate(ω, dt)` (exponential-map orientation integration for the rigid-body solver). Used for all scene rotation; glTF round-trips through this. Since tier 0 it also carries the ROBUST NORM pair: `magnitude()`/`normalise()` compute `sqrt(dot(q, q))` first and fall back to a scaled form whenever that sum is not finite and normal (zero, subnormal, infinite or NaN), which keeps ordinary rotations bit-identical while covering the ranges where the naive form has no accurate answer. `normalise` divides once per component wherever the length is representable — dividing twice costs an ulp each, which measurably changed solver behaviour — and its three answers are load-bearing: non-finite → a NaN quaternion (visibly invalid, NEVER laundered into the identity), below `float_normalise_cutoff` → the identity, otherwise normalised. `operator==` is exact component equality, so `q` and `-q` differ although they are the same rotation: the rotation-aware comparison is phase 2's job, not this type's. | | `math/mat3.hpp` | Column-major 3×3 (mirrors `Mat4`'s `[row,col]`): `fromQuaternion`, `diagonal`, `transpose`, `Mat3·Mat3` / `Mat3·Vec3`. Holds the world inverse inertia `R·diag(invI)·Rᵀ` in the physics solver. | | `math/mat4.hpp` | **Column-major.** Translation/rotation/scale/perspective/look-at. Everything downstream trusts this — verify the multiplication and handedness conventions. | | `math/view_basis.hpp` | Shared right/up construction that stays finite for zero-length or vertical look dirs. Used by view, skybox, shadow-fit, sort-depth. | diff --git a/docs/roadmap.md b/docs/roadmap.md index 74e7cee1..dce0cc78 100644 --- a/docs/roadmap.md +++ b/docs/roadmap.md @@ -22,7 +22,7 @@ the open items so they can't fork: | Doc | What it is | Status of its items | |---|---|---| -| [`codereview.md`](codereview.md) | Rolling **tiered static review**, following the [`review-order.md`](review-order.md) tiers (Tier 0 math, 18 Jul 2026; Tier 1 handles/limits/tunables, 19 Jul 2026). Further tiers expected. | **All open** — arc 3 below | +| [`codereview.md`](codereview.md) | Rolling **tiered static review**, following the [`review-order.md`](review-order.md) tiers (Tier 0 math, 18 Jul 2026; Tier 1 handles/limits/tunables, 19 Jul 2026). Further tiers expected. | Tier 0's correctness foundation is cleared (findings 1, 2, 4, 5 — see its § Phase 1 resolution); everything else open — arc 3 below | | [`architecturalreview.md`](architecturalreview.md) | One-shot **architectural review** (25 Jul 2026) of rendering, shadows/AA, physics, simplifier/VDPM. Audited 26 Jul so every finding now maps to a §6 row or an explicit "informational" tag. **Retire it once reviewed** — arc 2 below is self-contained. | 8 of 19 landed; the rest is arc 2 | | [`shadowplans.md`](shadowplans.md) | The **shadow-LOD improvement plan** (SH-01…SH-09) spun out of the architectural review's §2. | Milestones 0–2 landed (SH-01…SH-03, SH-05…SH-07, SH-04's deformation half); what remains is the follow-ups those left, plus evidence-gated SH-08/SH-09 — arc 1 | @@ -124,13 +124,18 @@ Handled as a unit the way CR-01…26 was, one branch per phase. Findings map to **Tier 0 — math & value types.** 3 high (`Mat3::inverse()` rejects valid small transforms; `approxEqual()` accepts NaNs as equal; rotation quaternions don't enforce their invariant), 5 medium (non-robust norms, "bitwise equality" isn't bitwise, affine/projective mixed, hidden projection -conventions, duplicated conversion authority) + a standardisation list. Sequenced by the doc: -1. **Correctness foundation** — NaN/tiny-matrix regression tests, shared scalar comparison, robust - scaled norms, scale-aware `Mat3::tryInverse()` + caller migration, fix/remove the bitwise API. -2. **Rotation redesign** — `UnitQuaternion`/`Rotation3`, one quaternion→matrix authority, migrate - transform/animation/render/physics users. *(Touches physics orientation ⇒ expect a determinism - golden re-baseline on BOTH platforms — see CLAUDE.md § Testing.)* -3. **Transform & API redesign** — `Affine3` + direct TRS, split affine point/vector/normal from +conventions, duplicated conversion authority) + a standardisation list. **Open here: the rotation +invariant, affine/projective separation, projection conventions and the conversion authority** — the +correctness foundation is done and its detail lives in [`codereview.md`](codereview.md), not in this +index. Sequenced by the doc: +1. **Rotation redesign** — `UnitQuaternion`/`Rotation3`, one quaternion→matrix authority, migrate + transform/animation/render/physics users. *(Touches physics orientation, so the goldens MAY move + — and "may" is the word. Do not pre-authorise a re-baseline: if a hash moves, identify the first + changed operation, confirm `ReplayIsBitIdentical` still holds, check the settle bound in + `Demos.Sleep`, and only then decide whether the new trajectory follows necessarily from enforcing + the rotation invariant. Phase 1 expected its goldens to move too, and the change that moved them + turned out to be a regression.)* +2. **Transform & API redesign** — `Affine3` + direct TRS, split affine point/vector/normal from projective, explicit projection conventions, standardise the access/operator surface. **Tier 1 — handles, limits, tunables.** 3 high open (texture generations not enforced on diff --git a/include/fire_engine/math/quaternion.hpp b/include/fire_engine/math/quaternion.hpp index 5389741d..00d128c4 100644 --- a/include/fire_engine/math/quaternion.hpp +++ b/include/fire_engine/math/quaternion.hpp @@ -148,9 +148,11 @@ class Quaternion } [[nodiscard]] - // Fast path first — see `VecBase::magnitude`. While the sum of squares is finite and positive - // this is the arithmetic the engine always did, bit for bit; the scaled form below runs only - // when that sum came back zero, infinite or NaN, which is precisely when it had no answer. + // Fast path first — see `VecBase::magnitude`. While the sum of squares is finite and NORMAL + // this is the arithmetic the engine always did, bit for bit; the scaled form below runs + // whenever it is not — zero, subnormal, infinite or NaN — which is precisely when it had no + // accurate answer. Subnormal is in that list because such a sum is finite and positive and has + // already lost most of its precision, which is the trap a `> 0` guard falls into. float magnitude() const noexcept { const float sumOfSquares = magnitudeSquared(); @@ -434,7 +436,8 @@ class Quaternion } private: - // THE ROBUST PATH, reached only when the sum of squares was zero, infinite or NaN. + // THE ROBUST PATH, reached whenever the sum of squares was not finite and NORMAL — zero, + // subnormal, infinite or NaN. [[nodiscard]] float scaledMagnitude() const noexcept { const float components[4]{x_, y_, z_, w_}; diff --git a/include/fire_engine/math/vec_base.hpp b/include/fire_engine/math/vec_base.hpp index d8d33892..e9759920 100644 --- a/include/fire_engine/math/vec_base.hpp +++ b/include/fire_engine/math/vec_base.hpp @@ -155,7 +155,7 @@ class VecBase } // Normalised through the same fast path, for the same reasons: while the sum of squares is - // finite and positive, this is exactly the division the engine did before — one rounding per + // finite and NORMAL, this is exactly the division the engine did before — one rounding per // component, bit-identical results. Rounding twice instead (dividing by the largest component // and then by a scaled norm) costs an ulp per component, which sounds like nothing and delayed // a settling box stack from step 169 to 425 on macOS and 1309 on Linux against a 600-step @@ -247,8 +247,10 @@ class VecBase } protected: - // THE ROBUST PATH, reached only when the sum of squares was zero, infinite or NaN — i.e. when - // the naive computation had no answer to give. Kept out of line from the hot path above. + // THE ROBUST PATH, reached whenever the sum of squares was not finite and NORMAL — zero, + // subnormal, infinite or NaN — i.e. whenever the naive computation had no accurate answer to + // give. Subnormal belongs in that list: such a sum is finite and positive and has already lost + // most of its precision. Kept out of line from the hot path above. [[nodiscard]] float scaledMagnitude() const noexcept { float largest = 0.0f; diff --git a/tests/physics/test_demos.cpp b/tests/physics/test_demos.cpp index 5f9f34e4..35844b6e 100644 --- a/tests/physics/test_demos.cpp +++ b/tests/physics/test_demos.cpp @@ -9,6 +9,7 @@ // The authored numbers here are the shared source of truth with // assets/physics_demos/generate.py — keep the two in sync when a demo changes. +#include #include #include #include @@ -433,9 +434,49 @@ TEST_CASE("Demos.Sleep.StackSleepsThenWakesOnImpact", "[Demos][slow]") } CHECK(anyAwake); - // Clean end state: the disturbance damps out and everything — stack and striker — - // comes to rest and sleeps on the floor. - step(world, 450); // -> step 600 + // Clean end state: the disturbance damps out and everything — stack and striker — comes to rest + // and sleeps on the floor. + // + // HOW LONG that takes is asserted, not just that it happens by some generous deadline. An + // endpoint check at step 600 passed on macOS while a change had tripled the settle time, and + // only failed on Linux because its trajectory was slower still: the tier-0 robust-norm work + // normalised through two divisions instead of one, costing an ulp per component, and friction + // directions that no longer cancel kept feeding the contact solver. Settle went 169 -> 425 here + // and 167 -> 1309 on Linux. A deadline generous enough never to flake is also generous enough + // to hide that, so the step count itself is the evidence. + // + // The bound is measured, not chosen for comfort: healthy runs settle at ABSOLUTE step 167 + // (macOS/arm64) and 172 (Linux/x86_64), so 260 leaves better than 50% headroom over the slower + // platform while sitting far below the 425 that a real regression produced. Widen it only with + // a measurement that says the engine legitimately got slower, and never to make a red test + // green. + // + // COUNTED IN ABSOLUTE STEPS, deliberately. Every number above and in the commit history is an + // absolute step index, and a loop counting from the impact window would silently mean something + // else: 260 iterations starting after step 150 is a bound of 410, which the 300-step regression + // this test exists to catch would pass comfortably. + constexpr int kImpactWindowEnd = 150; // the 90 + 60 stepped above + constexpr int kSettleByStep = 260; + int settledAtStep = -1; + for (int absoluteStep = kImpactWindowEnd + 1; absoluteStep <= kSettleByStep; ++absoluteStep) + { + step(world, 1); + const bool islandAsleep = + std::ranges::all_of(stack, [&](PhysicsBodyHandle h) { return world.sleeping(h); }) && + world.sleeping(striker); + if (islandAsleep) + { + settledAtStep = absoluteStep; + break; + } + } + INFO("the whole island slept at absolute step " << settledAtStep << " (budget " << kSettleByStep + << ")"); + CHECK(settledAtStep > 0); + + // And it STAYS asleep — a body that wakes itself again has not come to rest, it is oscillating + // around the threshold. + step(world, 120); for (const PhysicsBodyHandle h : stack) { CHECK(world.sleeping(h));