Unit quaternions have several advantages over 3x3 orthogonal matrices for representing 3d rotations on a computer.
However, one thing that has been disappointing me about the unit quaternion
representation is that axis-aligned 90 degree rotations
aren't exactly representable. For example, a 90-degree rotation around the z axis, taking the +x axis to the +y axis, is represented as [w=sqrt(1/2), x=0, y=0, z=sqrt(1/2)].
Surprising/unpleasant consequences include:
- applying a floating-point-quaternion-represented axis-aligned 90 degree rotation to a vector v often doesn't rotate v by exactly 90 degrees
- applying a floating-point-quaternion-represented axis-aligned 90 degree rotation to a vector v four times often doesn't yield exactly v
- squaring a floating-point-quaternion representing a 90 degree rotation around a coordinate axis doesn't exactly yield the (exactly representable) 180 degree rotation around that coordinate axis, and raising it to the eighth power doesn't yield the identity quaternion.
Because of this unfortunate lossiness of the quaternion representation on "nice" rotations, I still sometimes choose 3x3 matrices for applications in which I'd like axis-aligned 90 degree rotations, and compositions of them, to be exact and floating-point-roundoff-error-free. But the matrix representation isn't ideal either, since it loses the sometimes-needed double covering property (i.e. quaternions distinguish between the identity and a 360-degree rotation, but 3x3 rotation matrices don't) as well as other familiar desirable numerical properties of the quaternion representation, such as lack of need for re-orthogonalization.
My question: is there a computer representation of unit quaternions that does not suffer this imprecision, and also doesn't lose the double covering property?
One solution I can think of is to represent each of the 4 elements of the quaternion as a pair of machine-representable floating-point numbers [a,b], meaning a + b √2. So the representation of a quaternion would consist of eight floating-point numbers. I think that works, but it seems rather heavyweight; e.g. when computing the product of a long sequence of quaternions, each multiplication in the simple quaternion calculation would turn into 4 floating-point multiplications and 2 floating-point additions, and each addition would turn into 2 floating-point additions. From the point of view of trying to write a general-purpose library implementation, all that extra computation and storage seems pointless as soon as there's a factor that's not one of these "nice" rotations.
Another possible solution would be to represent each quaternion q=w+xi+yj+zk as a 4-tuple [sign(w)*w2, sign(x)*x2, sign(y)*y2, sign(z)*z2]. This representation is concise and has the desired non-lossiness for the subgroup of interest, but I don't know how to multiply together two quaternions in this representation.
Yet another possible approach would be to store the quaternion q2 instead of the usual q. This seems promising at first, but, again, I don't know how to non-lossily multiply two of these representions together on the computer, and furthermore the double-cover property is evidently lost.