Why are quaternions used for rotations?

Question

I'm a physicist, and have been learning some programming, and have come across a lot of people using quaternions for rotations instead of writing things in matrix/vector form.

In physics, there are very good reasons we don't use quaternions (despite the bizarre story that's occasionally told about Hamilton/Gibbs/etc). Physics requires that our descriptions have good analytic behavior (this has a precisely defined meaning, but in some rather technical ways that go far beyond what's taught in normal intro classes, so I won't go into any detail). It turns out that quaternions don't have this nice behavior, and so they aren't useful, and vectors/matrices do, so we use them.

However, restricted to rigid rotations and descriptions that do not use any analytic structures, 3D rotations can be equivalently described either way (or a few other ways).

Generally, we just want a mapping of a point X = (x, y, z) to a new point X' = (x', y', z') subject to the constraint that X² = X'². And there are lots of things that do this.

The naive way is to just draw the triangles this defines and use trigonometry, or use the isomorphism between a point (x, y, z) and a vector (x, y, z) and the function f(X) = X' and a matrix MX = X', or using quaternions, or projecting out components of the old vector along the new one using some other method (x, y, z)^T.(a,b,c) (x',y',z'), etc.

From a math point of view, these descriptions are all equivalent in this setting (as a theorem). They all have the same number of degrees of freedom, the same number of constraints, etc.

So why do quaternions seem to preferred over vectors?

The usual reasons I see are no gimbal lock, or numerical issues.

The no gimbal lock argument seems odd, since this is only a problem of euler angles. It is also only a coordinate problem (just like the singularity at r=0 in polar coordinates (the Jacobian looses rank)), which means it is only a local problem, and can be resolved by switching coordinates, rotating out of the degeneracy, or using two overlapping coordinate systems.

I'm less sure about numerical issues, since I don't know in detail how both of these (and any alternatives) would be implemented. I've read that re-normalizing a quaternion is easier than doing that for a rotation matrix, but this is only true for a general matrix; a rotation has additional constraints that trivializes this (which are built into the definition of quaternions) (In fact, this has to be true since they have the same number of degrees of freedom).

So what is the reason for the use of quaternions over vectors or other alternatives?