Things simplify considerably if you're willing to use quaternions.
In this case you have to simply store a "direction" unit vector(D) and an "upside" unit vector(U). Optionally you can store the "right" unit vector that points to the right direction (R) or you can compute doing the cross product (D x U).
When you want to change your pitch or bearing or orientation you simply have to construct a quaternion using the right vector and multiply the other two for that quaternion.
With this strategy you can store a "direction" vector that even points where you are going and rotate without the problem you may encounter if you touch some nasty rotation pole.
Quaternions are 4-tuple of real numbers that represent a real part and the coefficients for tree orthogonal imaginary parts i, j, k. A quaterion may be write as (s,v) where s is the real part and v is the imaginary 3D vector.
A point in 3D (p) may be represented as the quaternion (0,p) and the rotation of θ rads about the r unit vector is: q = (s,v) where s = cos(θ/2) and v is r·sin(θ/2).
now p* = q · p · q' is the rotated point where q' is the inverse of q.
The inversion of a unit quaternion equals its conjugate so q' = (s, -v)
The quaternion multiplication shows a little complexity: while the a quaternion can be written as q = w + x i +y j +z k thus p · q is a simple polinomial multiplication, there are the relations so:
we can say that a multiplication between quaternion can be written as:
this lead us to the conclusion that if q1 is the vector Q = [W1, X1, Y1, Z1] then
so q1 · q2 = Q · A in term of matrix multiplications.
You can see that a rotation using a quaternion involves 32 multiplications and 24 sums Vs the 9 multiplications and 6 sums for the rotation using the standard matrix rotations.
This difference is due to the fact that a fourth dimension is used to avoid the pole problem that often makes the matrix rotation unviable: the euclideans matrix rotations is not faster, it simply compute less things and you pay this if you pass through a pole.