From f5781813519883fd6c299a2ee31523942ead0a61 Mon Sep 17 00:00:00 2001
From: Christopher Durham <cad97@cad97.com>
Date: Fri, 12 Jan 2024 22:54:09 -0500
Subject: [PATCH] Fix glm::is_normalized epsilon

The existing comparison bound of $\epsilon^2$ is improperly scaled for
testing an epsilon of the squared vector magnitude. Let $\epsilon$ be
our specified epsilon and $\delta$ be the permissible delta of the
squared magnitude. Thus, for a nearly-normalized vector, we have

$$\begin{align}
\sqrt{1 + \delta} &=  1 + \epsilon        \\
          \delta  &= (1 + \epsilon)^2 - 1 \\
          \delta  &=  \epsilon^2 + 2\epsilon
\text{ .}\end{align}$$

Since we only care about small epsilon, we can assume
that $\epsilon^2$ is small and just use $\delta = 2\epsilon$. And in
fact, [this is the bound used by GLM][GLM#isNormalized] (MIT license)
... except they're using `length` and not `length2` for some reason.

[GLM#isNormalized]: https://github.com/g-truc/glm/blob/b06b775c1c80af51a1183c0e167f9de3b2351a79/glm/gtx/vector_query.inl#L102

If we stick an epsilon of `1.0e-6` into the current implementation,
this gives us a computed delta of `1.0e-12`: smaller than the `f32`
machine epsilon, and thus no different than direct float comparison
without epsilon. This also gives an effetive epsilon of `5.0e-13`;
*much* less than the intended `1.0e-6` of intended permitted slack!
By doing a bit more algebra, we can find the effective epsilon is
$\sqrt{\texttt{epsilon}^2 + 1} - 1$. This patch makes the effective
epsilon $\sqrt{2\times\texttt{epsilon} + 1} - 1$ which still isn't
*perfect*, but it's effectively linear in the domain we care about,
only really making a practical difference above an epsilon of 10%.

TL;DR: the existing `is_normalized` considers a vector normalized if
the squared magnitude is within `epsilon*epsilon` of `1`. This is wrong
and it should be testing if it's within `2*epsilon`. This PR fixes it.

For absence of doubt, a comparison epsilon of $\texttt{epsilon}^2$ is
correct when comparing squared magnitude against zero, such as when
testing if a displacement vector is nearly zero.
---
 nalgebra-glm/src/gtx/vector_query.rs | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/nalgebra-glm/src/gtx/vector_query.rs b/nalgebra-glm/src/gtx/vector_query.rs
index f3d0a9ea..a0b9f621 100644
--- a/nalgebra-glm/src/gtx/vector_query.rs
+++ b/nalgebra-glm/src/gtx/vector_query.rs
@@ -41,7 +41,10 @@ pub fn is_comp_null<T: Number, const D: usize>(v: &TVec<T, D>, epsilon: T) -> TV
 
 /// Returns `true` if `v` has a magnitude of 1 (up to an epsilon).
 pub fn is_normalized<T: RealNumber, const D: usize>(v: &TVec<T, D>, epsilon: T) -> bool {
-    abs_diff_eq!(v.norm_squared(), T::one(), epsilon = epsilon * epsilon)
+    // sqrt(1 + epsilon_{norm²} = 1 + epsilon_{norm}
+    // ==> epsilon_{norm²} = epsilon_{norm}² + 2*epsilon_{norm}
+    // For small epsilon, epsilon² is basically zero, so use 2*epsilon.
+    abs_diff_eq!(v.norm_squared(), T::one(), epsilon = epsilon + epsilon)
 }
 
 /// Returns `true` if `v` is zero (up to an epsilon).