and how to actually use them.

incomplete, under construction

Rotations

In most cases in computer graphics when you use quaternion, or axis and angle, or 3x3 matrix, you actually want to deal with 3D rotations. In a good design, interface (rotation) should be abstracted from internal representation (quaternion or 3x3 matrix).
Let's think what operations that are defined on rotations :

• Rotations can be concatenated. That is, if you rotate by A and then by B, you can find C that rotation by C gives same effect.
• For every rotation there is inverse rotation, so that rotation followed by inverse is no rotation.
• Comparison operator, 2 rotations may be equal(close within specific epsilon) or different.
• Applying rotation to vectors, or converting to matrix.
• Sending rotation into your API of choice (OpenGL or DirectX).
• For two rotations A and B and time t there is angular velocity W that turns from A to B in time t
• Rotation can be mirrored on some axis.

and so on and so forth.

Concatenation:

It is reasonable to use operator*(Rotation, Rotation) for that. (Why not '+'? Because that'll screw up your algebra a lot. You have + already defined on vectors.)
So, (B*A) defines rotation by A followed by rotation by B (or, reference frame A as child of reference frame B. More on why this order is preferred below).
There you can see that it's like same old multiplication of matrices or quaternions.
So you can ask "what's the point in defining same thing again". What is very important to understand is that multiplication of matrices or quaternions do more than concatenation of rotations. For example, it can be used for scaling, or shearing.
Computational inprecision, over time, leads to de-orthonormalization of matrix, or to non-unit-length quaternion, which not only rotates but also distorts or scales the object.
You can add normalization of quaternions into this operator*(Rotation, Rotation), to ensure stability. Or you can add debugmode assertion of 1-epsilon>= |Q| <=1+epsilon, to simplify debugging of computational instability problems.

Inverse:

There things get much more interesting.
You dont need to really compute inverse of quaternion or 3x3 matrix to inverse rotation!
You only need transpose of 3x3 matrix, or negation of real part of quaternion (it's actually -conjugate(Q). Remeber that -Q stores same rotation as Q, so -conjugate(Q) and conjugate(Q) define same rotation).
While with matrices and quaternions your code would need to use different functions, transpose and conjugate, Rotation class could neatly hide and abstract those details as Inverse.
Note that it really does hide implementation details from rest of code (unlike accessors like get_X and set_X for instance which hide nothing).

Euclidean transform (almost)

That is all nice but usually you want something more than rotations. You want rotation and translation.
That is euclidean transform, sans the reflections. I called it CoordSys in my math library.
How do we do that with our Rotation class? Easily, store rotation and vector for translation.
There we do have a slight problem. What is multiplication and what is inverse for our euclidean transform class?
The transform is defined as:

v'=T+R*v

where T is translation vector and R is rotation.
Naturally, concatenation of transform T1 R1 followed by T2 R2 is

v'=T2+R2*(T1+R1*v) = T2 + R2*T1 + R2*R1*v

For composite rotation we have

T = T2+R2*T1
R = R2*R1

That is, we can write our multiplication operator as something like

template
CoordSys operator*(const CoordSys &c1, const CoordSys &c2)
{
    CoordSys result;
    result.orientation=
#ifndef dont_normalize
        Normalized
#endif
        (c1.orientation*c2.orientation);
    result.position=ToMatrix(c1.orientation)*c2.position+c1.position;
    return result;
};

Inverse:

The transform, again, is

v'=T+R*v

using some simple algebra on it to solve for v:

v'-T=R*v
inverse(R)*(v' - T) = inverse(R)*R*v = v
v = -inverse(R)*T +inverse(R)*(v')

and in code it is similarly trivial:

template
CoordSys Inverse(const CoordSys &c1)
{
    CoordSys result;
    result.orientation=Inverse(c1.orientation);
    result.position= -(result.orientation*c1.position);
    return result;
};

Now, you can use CoordSys in place of 4x4 matrix for your transformations, using quaternion for the rotation. It takes up over 2x less space than 4x4 matrix, it does not suffer from as many numerical problems (it can not store unwanted transformations such as shear), and quaternion rotations are easier to interpolate than matrices.

incomplete, under construction

Todo: clean up my implementation and publish (under BSD license).

---

References:

MathWorld: Quaternion.
EuclideanSpace.com: Maths - Combined Rotation and Translation (great site, BTW.)
http://gandalf-library.sourceforge.net/ some library that includes Euclidean transformations as well.
Geometric Transformations