In my initial code I had player movement applied directly as offsets to the position and orientation data.
My idea was to try and change that, so that player movement and rotation was instead applied once to the player's velocity or angular velocity vectors.
With simple translation this is easy. However, when the player rotates, the translation vector needs to change direction. This means subtracting the old vector and adding the new one. To complicate things further, the same is true of the 'up' vector if the player is pitching downwards.
So before each physics update the old vectors are subtracted, a new orientation matrix is calculated for the next frame, and the velocity vectors are re-applied in the proper directions.
This worked for either translation or rotation alone, but when moving and rotating simultaneously I ended up with residual velocity (around 0.02-0.07 units/sec, compared to a movement speed of around 1.6-8.0 units/sec). It's small enough that it seems like it could possibly be a rounding error but also not so small that there could be some small error in my algorithm, or else the method described above is just not possible.
I wanted just one set of velocity and angular velocity vectors, so that during the physics update step all collision could be handled in one step. Otherwise I will have to make a second collision step if I simulate player movement separately. This seems undesirable since my player is just an instance of a mobile game entity which I will be adding later, so all these objects would add more overhead to the movement step.