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/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/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 4a90948e..a962dadd 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 @@ -957,6 +961,38 @@ 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 + 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 + (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 b9fcb4c3..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). `operator==` is **strict bit equality** — use `approxEqual` for tolerance. 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/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. 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/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/mat3.hpp b/include/fire_engine/math/mat3.hpp index f38ec2ad..8d49799a 100644 --- a/include/fire_engine/math/mat3.hpp +++ b/include/fire_engine/math/mat3.hpp @@ -1,11 +1,25 @@ #pragma once +#include +#include +#include + #include +#include #include 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. @@ -101,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_) + { + if (!std::isfinite(value)) + { + return std::nullopt; + } + const double magnitude = + value < 0.0f ? -static_cast(value) : static_cast(value); + scale = magnitude > scale ? magnitude : scale; + } + if (scale == 0.0) { - return Mat3{}; + return std::nullopt; // the zero matrix: singular, and the one case scaling cannot help } - const float s = 1.0f / det; + + 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; } @@ -213,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 { @@ -227,13 +331,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 +350,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..ddce1482 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 @@ -94,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 { @@ -109,19 +113,18 @@ class Mat4 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 bitwiseEqual(const Mat4& rhs) const noexcept - { - return *this == rhs; - } - - [[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 +132,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..00d128c4 100644 --- a/include/fire_engine/math/quaternion.hpp +++ b/include/fire_engine/math/quaternion.hpp @@ -1,9 +1,11 @@ #pragma once #include +#include #include #include +#include #include namespace fire_engine @@ -85,31 +87,46 @@ 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_; } + // 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 bitwiseEqual(const Quaternion& rhs) const noexcept + constexpr bool approxEqual(const Quaternion& rhs, float eps, float relativeEps) const noexcept { - return *this == rhs; + 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 = float_epsilon) const noexcept + constexpr bool approxEqual(const Quaternion& rhs, float eps) 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; + // 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]] @@ -131,20 +148,51 @@ class Quaternion } [[nodiscard]] + // 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 { - 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 @@ -388,6 +436,82 @@ class Quaternion } private: + // 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_}; + 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/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..e9759920 100644 --- a/include/fire_engine/math/vec_base.hpp +++ b/include/fire_engine/math/vec_base.hpp @@ -2,8 +2,10 @@ #include #include +#include #include +#include namespace fire_engine { @@ -97,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; @@ -114,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 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 + // 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) + 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)) { - return Derived{}; - } - - Derived result{v}; - for (std::size_t i = 0; i < N; ++i) - { - 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 @@ -137,8 +196,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 { @@ -152,24 +214,131 @@ class VecBase 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 Derived& rhs, float eps, float relativeEps) const noexcept + { + for (std::size_t i = 0; i < N; ++i) + { + if (!almostEqual(data_[i], rhs.data_[i], eps, relativeEps)) + { + return false; + } + } + return true; + } + [[nodiscard]] - constexpr bool bitwiseEqual(const Derived& rhs) const noexcept + constexpr bool approxEqual(const Derived& rhs, float eps) const noexcept { - return self() == rhs; + // ABSOLUTE ONLY — a stated tolerance is the whole answer. + return approxEqual(rhs, eps, 0.0f); } [[nodiscard]] - constexpr bool approxEqual(const Derived& rhs, float eps = float_epsilon) const noexcept + constexpr bool approxEqual(const Derived& rhs) const noexcept { + return approxEqual(rhs, float_epsilon, float_relative_epsilon); + } + +protected: + // 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; + bool anyInfinite = false; for (std::size_t i = 0; i < N; ++i) { - const float diff = data_[i] - rhs.data_[i]; - if (diff > eps || diff < -eps) + if (std::isnan(data_[i])) { - return false; + 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; } - return true; + 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: 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()); } 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..0c170d54 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]") @@ -448,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_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)); +} diff --git a/tests/math/test_vec3.cpp b/tests/math/test_vec3.cpp index 526cc812..00d5c50b 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]") @@ -809,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{}); +} 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));