Not sure it was your question, but ...

Stability of the rotation axis :

If you multiply two quaternions around the same axis (N), you can see the multiplication formula geometrically. It probably shows why the resulting axis is well preserved numerically.

q1 * q2 = { S1^S2 + c1*S1 + c2*S1, c1*c2 - S1*S2 }
(S1=s1*N and S2=s2*N)

Since S1//S2//N, S1^S2 = 0 (1)
c1.S1 + c2.S2 is also // N. (2)

(1) In practice I suppose s1 and s2 are usually small => second order epsilon of small angles. The 'N'^'N' (S1^S2/(s1*s2)) part is also very small, second order of the preceding numerical unstabilities. Thus I suppose it should become a true 0 in practice when added to c1*S1 + c2*S1 (due to 24 bits precision).

(2) The multiplications preserve directions a lot. I would be a bit harder to show. I feel it ;)

As a consequence if you apply many subsequent rotations around the same axis, intuitively, it should keep very stable numerically.

@Eelco

Since I feel your issue is linked to the previous discussions we had : the kinetics and collision prediction. Your concern is interpolating body frames (ct angular velocities) and retrieving vertex positions. Dynamics do not interfere at the resolution of the small time intervals we consider, (except when collisions or contacts occur). Here is a proposal to limitate the tradeoffs of the previously mentionned solutions :

1) - any angular velocity
2) - numerical precision/stability
3) - code speed

Your algo starts with {Axis, W} (AW) structures as inputs. This solves 1). When a velocity is too high you know that our quadratic approximations won't work, so you need to slice the interval anyway.

For the manageable intervals or interval divisions, your AW can then be converted to a quaternion. Eventully with powers of 1/2, you have not lost any precision up to here. So 2) is OK.

For 3) you can check w (scaled to the time interval) before converting. You can avoid inverse trigonometrics for small angular speeds (w = max angle), since sinus(w/2) = w/2 + o(w^3). You just need a NR sqrt for cos (1-w*w/4 will do).

Also for such cases, if required, lerp + sqrt ( S'=t*S and c=sqrt(1-(t*t*s*s)) ) will be as precise as slerp.