Math correctness foundation - #147
Merged
Merged
Conversation
Every approximate comparison in the engine — VecBase, Mat3, Mat4,
Quaternion — was the same four lines:
const float diff = a - b;
if (diff > eps || diff < -eps) { return false; }
Both comparisons are false when diff is NaN, so a NaN compared
approximately equal to anything, itself included. Two equal infinities
subtract to NaN and passed the same way, which reads as correct until
you notice that opposite infinities do too. A test asserting a transform
had stayed finite would pass on a transform that was entirely NaN:
precisely the failures tests exist to catch, hidden by the thing
catching them.
math/scalar.hpp is now the only place that decides whether two numbers
are close. It asks `a == b` first, so identical values are equal
whatever they are — including equal infinities, which the old form got
right by accident rather than by rule. Any remaining non-finite operand
is unequal. Finite values compare against the larger of an absolute and
a relative tolerance, in double intermediates: `a - b` in float
overflows for values near FLT_MAX of opposite sign, so the difference
itself could be non-finite for two perfectly finite inputs and the
comparison would answer about a number neither caller mentioned.
THREE OVERLOADS, because a tolerance a caller spells out must mean what
it says. No tolerance takes both defaults; one tolerance is ABSOLUTE
ONLY, exactly the historical comparison; two are both terms, stated. A
single function with a defaulted relative term would have made
approxEqual(rhs, 1e-9f) admit a 1e-7 difference at magnitude 1.0 — an
explicit tolerance silently overridden by an implicit one. The existing
strictness tests are unchanged, which is the evidence that the overloads
did their job.
TOLERANCES ARE VALIDATED, not reinterpreted. A negative, NaN or infinite
tolerance returns false, and the check runs before the `a == b` shortcut
so equal operands cannot hide it. Treating a negative as "a term that
cannot be satisfied" is mathematically defensible and diagnostically
useless: it converts a broken configuration into a different valid
policy and lets the run continue looking healthy.
float_relative_epsilon is its own constant. float_epsilon was already
serving as absolute tolerance and degeneracy threshold; a third job
would have made it a number nobody could change safely.
Mat3::inverse() failed by returning the zero matrix — both a legitimate value and an error report, so a caller could not tell "singular" from "the answer is zero" without asking again. Its threshold was ABSOLUTE (|det| <= 1e-12), which rejects transforms that are perfectly well conditioned merely for being small: a uniform scale of 1e-5 has determinant 1e-15 and an exact inverse of 1e5, and the old test called it singular. VDPM had already worked around that with its own scale-invariant predicate, then called the absolute-threshold inverse anyway — so `coneUsable` reported true while `cameraObj` was a zero matrix times a camera offset. Two invertibility decisions, disagreeing, in adjacent lines. tryInverse() returns std::optional<Mat3> and is the only one. It normalises the matrix by its largest absolute component before taking the determinant, so the question asked is conditioning rather than size; it computes in double; it refuses non-finite input rather than producing a matrix of NaNs that every later operation spreads silently. Two things the optional means that a determinant test alone does not. An invalid tolerance is refused before the matrix is examined — `magnitude > tolerance` is TRUE for a zero determinant against a negative threshold, so a singular matrix would otherwise be accepted and then divided by its own zero. And every element is checked finite and within float's range in double BEFORE any conversion: a uniform scale of 1e-39 is perfectly conditioned and its inverse is 1e39, which float cannot hold, so converting would hand back an engaged optional full of infinities. The optional means "a representable inverse exists" or it means nothing. determinant() now returns double, and VDPM is why. A reflected instance with a tiny uniform scale has a usable inverse and a determinant that underflows float to -0.0f — which is not less than zero, so the facing sign says "winding preserved" for a transform that reverses it, and the cone culls the side that should be visible. The threshold for the underflow is lower than it looks: 1e-5 cubed is fine, 1e-16 cubed is not. VDPM's conditioning POLICY is reproduced exactly rather than inherited. tryInverse normalises by largest component while the old predicate used σ_max, and the two differ by (σ_max/maxComponent)³ ∈ [1, 27]; passing a constant would have relaxed the policy by up to 27× (or 27000× with the permissive default), and a blanket conservative constant would have REJECTED shapes the old test accepted. Scaling the tolerance by that exact ratio gives |det| > 1e-6·σ_max³ again — one decision, same answer, both sides of the boundary pinned by tests. The three physics inversions — the D block and Schur complement inside SpatialMatrix::inverse(), and invertDof's 3-DOF path — have a genuine invariant, so they go through a physics-local invertInvariant(): assert, and in release a logged reason and abort. An assert alone would leave *inverse dereferencing an empty optional in exactly the situation where the solver state is already corrupt. The helper stays in physics deliberately; in math/ it would recreate the unconditional inversion API this removes. Note for anyone reading the goldens: Determinism.GoldenHash builds a floor and three rigid boxes and never constructs an articulation, so it does not exercise these paths. The articulation suite is the relevant coverage.
bitwiseEqual() was `return self() == rhs` — ordinary floating-point comparison wearing a name that promised bit patterns. It was wrong in both directions it could be: -0.0f and +0.0f have different bits and compare EQUAL, while a NaN has bits identical to itself and compares UNEQUAL. It also duplicated operator== exactly, so it offered a choice between two spellings of one thing, one of which lied. Nothing outside math/ ever called it. Deleted rather than corrected: a method that means what the operator means earns its place only by being clearer, and this one was less clear. The comments claiming "strict bit-for-bit equality" are corrected on all four types — including Mat3, which made the claim without the method — and now state the two IEEE cases outright, plus what a caller would actually need if a determinism diagnostic ever wants real bit comparison (std::bit_cast, explicitly, at the point that needs it). The three tests were named BitwiseEqualMatchesOperator and asserted that the method agreed with the operator it forwarded to, which could not fail. They now assert the semantics the corrected comments claim: signed zeros compare equal, NaNs do not compare equal to themselves, and — for Quaternion — q and -q compare unequal despite being the same rotation, which is the seam the rotation redesign in phase 2 has to close.
sqrt(dot(v, v)) fails at both ends. Components above ~1.8e19 square to infinity, so a vector whose length is perfectly representable reports infinity and then normalises to zero; components below ~1e-22 square to zero, so a small but ordinary vector reports a length of zero and loses its direction entirely. Neither case announces itself as a failure — both return a plausible number. Both failures DO announce themselves in the sum, which comes back infinite, NaN, zero, or subnormal. So the naive sum is computed first and trusted whenever it is finite and NORMAL, and the scaled fallback runs only where there was no correct answer before it. That ordering is the whole design: - Ordinary vectors take the arithmetic the engine always used, bit for bit. Neither physics golden moves. A change about extreme values has no business perturbing a box stack, and this one does not. - The robust path is confined to inputs that were already wrong. Normal, not merely positive, and that distinction is the difference between covering the underflow region and covering only its endpoint. A sum that has gone subnormal is finite and nonzero and has already lost most of its precision: (3e-23, 3e-23, 0) sums to 2.8e-45, carries about two significant bits, and answered 5.29e-23 against a true 4.24e-23 — 24.8% high, from a fast path that believed itself safe because the sum looked usable. A finite subnormal sum does reach the scaled fallback in normalise() as well; it simply cannot produce a normalised result, because such a magnitude is necessarily below float_normalise_cutoff and the configured degenerate value comes back instead. normalise() divides ONCE per component wherever the length is representable. The first version of this commit divided twice — by the largest component, then by a scaled norm — which is an extra ulp per component and looked harmless. It delayed a settling box stack from step 169 to 425 on macOS and from 167 to 1309 on Linux, against a 600-step budget, so the sleep demo failed on Linux only and looked for all the world like a marginal test. It was not marginal: friction directions that no longer cancel keep feeding the contact solver. A solver notices an ulp; that is what a solver is. The two-step form survives only for the case that needs it — a vector of finite components whose norm exceeds float's range, where dividing by the (correctly infinite) length would answer zero and throw away an ordinary direction. Three answers from normalise, and the middle one is the one that matters: a non-finite component yields NaNs, VISIBLY invalid; a magnitude below float_normalise_cutoff yields the zero vector or the identity rotation, as it always has; anything else is normalised. Laundering a corrupt input into the degenerate answer would turn a NaN orientation into a confident "no rotation" and surface it three seconds later somewhere unrelated. float_normalise_cutoff is its own constant at the same numeric value float_epsilon had. One number was serving as absolute comparison tolerance and degeneracy threshold at once; the rename moves no behaviour and stops the next person changing one policy while meaning the other. Measured, 2M calls: Quaternion::normalise 2.39 -> 2.49 ms (parity), Vec3::magnitude 3.27 -> 5.59, Vec3::normalise 2.68 -> 4.07. The residual on Vec3 is the guard branch defeating vectorisation in a synthetic back-to-back loop; swapping isfinite for a bare comparison measured identically (5.67 vs 5.59), so the clearer form stayed. Engine-level physics workloads are unchanged.
Demos.Sleep checked that everything was asleep by step 600. That deadline was generous enough never to flake and therefore generous enough to hide a regression: when the tier-0 norm work added an ulp of error per component, the settle time went from step 169 to 425 on macOS and the test still passed, with 175 steps to spare. Only Linux failed, because its trajectory was slower still — 1309 — and it looked exactly like a marginal test on a slow platform. The obvious response, widening the budget, would have buried the defect permanently. The step count is now the evidence. The test records the first post-impact step at which the whole island — the three-box stack AND the striker — is asleep, and bounds it at 260. That number is measured, not chosen for comfort: healthy runs settle at 167 (macOS/arm64) and 172 (Linux/x86_64), so 260 leaves better than 50% headroom over the slower platform while sitting decisively below the 425 a real regression produced. The gap between those two figures is what makes the assertion worth having; a bound inside it can distinguish them. It then steps a further 120 and re-checks, which separates "crossed the sleep threshold once" from "came to rest". A body oscillating around the threshold satisfies the first and fails the second. Widen the bound only with a measurement showing the engine legitimately got slower, and never to make a red test green — the comment in the test says so, because that is exactly the pressure this assertion will apply. The count is in ABSOLUTE simulation steps, and that is not incidental. Every figure here and in the branch history — 167, 172, 425, 1309 — is an absolute step index, so a loop counting from the end of the impact window would mean something quite different: 260 iterations starting after step 150 is a bound of 410, which is above the regression this test exists to catch. It reported settling at 19 where every measurement says 169. The loop therefore runs from kImpactWindowEnd + 1 to kSettleByStep and reports the absolute step. Also closes tier-0 phase 1 in the docs: docs/codereview.md records the resolution against findings 1, 2, 4 and 5 as five commits (this one included, since it exists only because of how the norms regression was caught), and docs/roadmap.md drops phase 1 entirely — it indexes open work, and this phase is no longer open.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.