> they only work if the quaternion's magnitude is exactly one

That's why you always normalize the quaternion first, and the article seems to require the normalized form:

Q.54 How do I convert a quaternion to a rotation matrix?

Assuming that a quaternion has been created in the form:

Q = |X Y Z W|

At least, I would read |X Y Z W| as meaning normalized(X Y Z W)

I don't see this notation explicitly defined when they describe quaternion normalization (Q.52) though, so I agree this leaves much out. It's more a cheat sheet than learning material.

> You need to replace sqrt(1 - a*a) assumptions with actual components, and use atan2 instead of acos

I'm kind of rusty with this, but I think the reason we don't do that is that it's cheaper to normalize then convert rather than use the non-normalized conversion formula. Correct me if I'm wrong.

> That's why you always normalize the quaternion first

Again, that won't solve the problem of floating point components not perfectly adding up to one.

Constantly normalizing quats are unnecessary performance hits and worsens your accuracy with no benefit other than being marginally "easier for low-math programmers to reason about", when you could instead just work with non-unit quats homogeneously and divide the final result with the quat square magnitude just once in the entire pipeline, or not at all and instead just pass the square magnitude as your w factor that gets divided on the GPU anyways