Oh, you're right. Thanks! That was a good hint. I gotta go change some code...
Representing uniform scaling in a quaternion works nicely when the conversion as layed out here
is used (the matrix representation right before the sentence "since qw*qw + qx*qx+ qz*qz+ qy*qy = 1 this gives").
Apparently, one also does not need to renormalize the quaternion when applying any rotations. This is really nice. Applying uniform scaling to any quaternion is then just scalar-multiplying by sqrt(scaleFactor).
Yeap. Since rotation is P' = QPQ-1
(which eliminates any scale factors) or QPQ*
which will scale by magnitude of Q squared and these compose as expected.