After seeing this thread, I was trying to think of how I would pack a 3D rotation into a fixed number of bits, for instance to store animation data. I wrote a little piece of code with the result and I think it's neat enough to share.

The idea is to encode a unit-length quaternion in a manner analogous to cube mapping, but with one extra dimension and taking advantage of the property that a sign flip doesn't change which rotation is being represented.

One can think of cube mapping as consisting of an encoding of points in a sphere by first indicating which coordinate has largest absolute value and what sign it has (i.e., which of the 6 faces of the axis-aligned cube the point projects to) and then the remaining coordinates divided by the largest one (i.e., what point of the face the point projects to).

In our case, we don't need to encode the sign of the largest component, so we only need to use 2 bits to encode what the largest component is, and we can use the remaining bits to encode the other three components.

I think 32 bits is probably good enough for animation data in a game, and it's convenient that 30 is a multiple of 3, so it's easy to encode the other components. Actually, even if we didn't have that convenience, it wouldn't be a big deal to use a resolution that is not a power of 2, but some integer divisions would be involved in the unpacking code.

Here's the code, together with a main program that generates random rotations and measures how bad the dot product between the original and the packed and unpacked gets (the dot product seems to be > 0.999993, although I haven't made a theorem out of it):
#include <iostream>
#include <cstdlib>
#include <boost/math/quaternion.hpp>

typedef boost::math::quaternion<double> quaternion;

int double_to_int(double x) {
return static_cast<int>(std::floor(0.5 * (x + 1.0) * 1023.0 + 0.5));
}

double int_to_double(int x) {
return (x - 512) * (1.0 / 1023.0) * 2.0;
}

struct PackedQuaternion {
// 2 bits to indicate which component was largest
// 10 bits for each of the other components
unsigned u;

PackedQuaternion(quaternion q) {
int largest_index = 0;
double largest_component = q.R_component_1();
if (std::abs(q.R_component_2()) > std::abs(largest_component)) {
largest_index = 1;
largest_component = q.R_component_2();
}
if (std::abs(q.R_component_3()) > std::abs(largest_component)) {
largest_index = 2;
largest_component = q.R_component_3();
}
if (std::abs(q.R_component_4()) > std::abs(largest_component)) {
largest_index = 3;
largest_component = q.R_component_4();
}

q *= 1.0 / largest_component;

int a = double_to_int(q.R_component_1());
int b = double_to_int(q.R_component_2());
int c = double_to_int(q.R_component_3());
int d = double_to_int(q.R_component_4());

u = largest_index;
if (largest_index != 0)
u = (u << 10) + a;
if (largest_index != 1)
u = (u << 10) + b;
if (largest_index != 2)
u = (u << 10) + c;
if (largest_index != 3)
u = (u << 10) + d;
}

quaternion get() const {
int largest_index = u >> 30;
double x = int_to_double((u >> 20) & 1023);
double y = int_to_double((u >> 10) & 1023);
double z = int_to_double(u & 1023);

quaternion result;
switch (largest_index) {
case 0:
result = quaternion(1.0, x, y, z);
break;
case 1:
result = quaternion(x, 1.0, y, z);
break;
case 2:
result = quaternion(x, y, 1.0, z);
break;
case 3:
result = quaternion(x, y, z, 1.0);
break;
}

return result * (1.0 / abs(result));
}
};

double rand_U_0_1() {
return std::rand() / (RAND_MAX + 1.0);
}

quaternion random_rotation() {
quaternion result;
do {
result = quaternion(rand_U_0_1()*2.0-1.0, rand_U_0_1()*2.0-1.0, rand_U_0_1()*2.0-1.0, rand_U_0_1()*2.0-1.0);
} while (norm(result) > 1.0);
return result*(1.0/abs(result));
}

double dot_product(quaternion q, quaternion p) {
return q.R_component_1() * p.R_component_1() +
q.R_component_2() * p.R_component_2() +
q.R_component_3() * p.R_component_3() +
q.R_component_4() * p.R_component_4();
}

int main() {
double worst_dot_product = 1.0;
for (int i=0; i<1000000000; ++i) {
quaternion q = random_rotation();
PackedQuaternion pq(q);
quaternion p = pq.get();
if (dot_product(p,q) < 0)
p *= -1.0;
if (dot_product(p,q) < worst_dot_product) {
worst_dot_product = dot_product(p,q);
std::cout << i << ' ' << q << ' ' << p << ' ' << worst_dot_product << '\n';
}
}
}


Any comments are welcome, and feel free to use the idea or the code if you find them useful.