The PR is failing the unit tests. Please, merge master to your branch and resolve the failure for the PR to be merged.

Thanks,
Christophe

So, I have investigated this a bit further.

I had gone a bit fast with the threshold value. I wrote a small program to iterate over all possible float32 angles from 0 to π / 2 to compare the loss of precision of the cos-method and the sin-method. It seems like the sin-method is worth it for angles roughly below 1 rad, and the cos-method for angles above this value.

This means the value w should be compared to cos(0.5) ≈ 0.878, not to 0.5. The unit test which was failing uses an angle of π/2, or a quaternion with w ≈ 0.707, and expects an exact value (the result, which lies between 1.0 and 2.0, is compared with an absolute strict difference of one epsilon).

Also, the sin-method is not slower than the cos-method on the [0.0, 1.0] interval, at least with a naive scalar implementation.

So, the new version passes the unit tests, and should essentially extract most of the available precision, without sacrificing much performance.

To visualize better what this PR does, this is the imprecision of acos(cos(x)):

We easily notice the loss of precision in [-1; 1]. The change in behavior looks sharp, but do note that the x-axis is given in logarithmic scale. For |x| > 1, ULP(x) ≥ ε; since |cos(x)| ≤ 1, it does not waste precision. However, for |x| < 1, ULP(x) < ε (and particularily for |x| ≪ 1, ULP(x) ≪ ε); in this case ½ ≤ |cos(x)| < 1, so ULP(cos(x)) = ε / 2. It means that, as x gets close to 0, the absolute precision of its float32 representation increases, while that of cos(x) remains constant. In layman's terms, more and more possible values for x gets mapped to the same value by cos.

Outside the ±1 boundary, you get most of the precision back, so the curve looks almost flat. Within it, you quickly gets >100 % error, so it plateaus at what looks like 62 bits of error (maybe 64-bit fixed-point representation, 1 bit of sign, and 1 bit because a broken clock is right twice a day?). Between roughly ±1e-19 all precision is lost (even with 64 bits of precision), and acos(cos(x)) always return exactly 0. However, the relative precision gets lower, as the real values gets closer and closer to the approximated 0.

Compare this to the curve for asin(sin(x)):

Because sin(x) ~ x around 0, we observe no large loss of precision around there. Outside ±1, we have the same behavior as acos(cos(x)).

Although asin(sin(x)) behaves much more nicely than acos(cos(x)), the sin-method of retrieving the angle of a quaternion also implies computing √(x²+y²+z²). The loss of precision of the full sin-method is shown below (note that it actually depends on the axis used):

Clearly, we are still losing quite a few bits of precision. So it's best to use the cos-method outside ±1, and the sin-method inside ±1.

Using hypot can improve things by replacing sqrt(x*x + y*y + z*z) with hypot(hypot(x, y), z). The result looks almost too good:

By manually iterating over the floats, it seems like the sin-hypot method indeed behaves much more nicely than the sin-sqrt method for values lower than roughly 2⁻⁶². Actually, it behaves well even on subnormal numbers, and only starts to break down below 2⁻¹²⁶. My guess is that this might be related to the 64 bits of the extended precision of the FPU. However, the asin-hypot method is several times slower than the asin-sqrt method where it matters. It might make sense to use the asin-hypot method below a set threshold, but the appropriate value for this threshold would most likely depend on the architecture.

Finally placated the gods of the CI 😅

I really wish all the PR would be of that level of quality !

Thanks for contributing !
Christophe

Thank you very much for your patience and your kinds words. GLM is really a great library, I was happy to be able to contribute!