Computer representation of quaternions that's exact for 90-degree rotations?

Question

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)*w², sign(x)*x², sign(y)*y², sign(z)*z²]. 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 q² 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.

Answer 1

You probably want to check the paper "Algorithms for manipulating quaternions in floating-point arithmetic" published in 2020 available online:

https://hal.archives-ouvertes.fr/hal-02470766/file/quaternions.pdf

Which shows how to perform exact computations and avoid unbounded numerical errors.

EDIT:

You can get rid of the square roots by using an unnormalized (i.e., non-unit) quaternion. Let me explain the idea.

Having two unit 3D-vectors X and Y, represented as pure quaternions, the quaternion Q rotating X to Y is

Q = X * (X + Y) / |X * (X + Y)|

The denominator, which is taking the norm, is the problem since it involves a square root.

You can see it by expanding the expresion as:

Q = (1 + X * Y) / sqrt(2 + 2*(X • Y))

If we replace X = i and Y = j to see the 90 degree rotation you get:

Q = (1 + k) / sqrt(2)

So, |1 + k| = sqrt(2). But we can actually use the unnormalized quaternion Q = 1 + k to perform rotations, all we need to do is to normalize the rotated result by the SQUARED norm of the quaternion.

For example, the squared norm of Q = 1 + k is |1 + k|^2 = 2 (and that is exact as you never took the square root) lets apply the unnormalized quaternion to the vector X = i:

= (1 + k) i (1 - k)

= (i + k * i - i * k - k * i * k)

= (i + 2 j - i)

= 2 j

To get the correct result we divide by squared norm.

I haven't tested but I believe you would get exact results by applying unnormalized quaternions to your vectors and then dividing the result by the squared norm.

The algorithm would be

1. Create the unnormalized quaternion Q = X * (X + Y)

2. Apply Q to your vectors as: v' = Q * v * ~Q

3. Normalize by squared norm: v'' = v' / |Q|^2