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 X2 = X'2. 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?
Gimbal lock is one reason, although as you say it is only a problem with Euler angles and is easily solvable. Euler angles are still used when memory is a concern as you only need to store 3 numbers.
For quaternions versus a 3x3 rotation matrix, the quaternion has the advantage in size (4 scalars vs. 9) and speed (quaternion multiplication is much faster than 3x3 matrix multiplication).
Note that all of these representations of rotations are used in practice. Euler angles use the least memory; matrices use more memory but don't suffer from Gimbal lock and have nice analytical properties; and quaternions strike a nice balance of both, being lightweight, but free from Gimbal lock.