As an alternative to use of quaternion math:
[This method requires ability to do a generalized coordinate rotation in 3D computationally].
There is a fairly simple way to convert from [ Euler axis, Euler angle ] representation into a 3x3 rotator matrix. Once the two rotators are in matrix form, you just bang them together into the composite transform R1 • R2.
To convert [ Euler axis, euler angle ] to 3x3 matrix:
Transform the standard x-axis direction [ 1 0 0 ] using the Euler spec.
a. First put direction [ 1 0 0 ] into rotated coordinates where the newZaxis takes on the value of the Euler axis direction.
b. Call a rotation of the Euler angle amount about the z-axis on the results from a.
c. Call the inverse rotation transform used in a. on the results from b.
You now have a direction vector giving the value of newXaxis (first column of your 3x3 rotator).
Repeat the same transform process on the other standard axes:
process [ 0 1 0 ] to obtain the newYaxis
process [ 0 0 1 ] to obtain the newZaxis
Finally compose your equivalent 3x3 rotator matrix R
R = [ newXaxis newYaxis newZaxis ]
and use matrix composition to combine any 2 two such rotators.
The final part of your questions appeals to how to convert a 3x3 rotator matrix into [ Euler axis, Euler angle ] representation.
I've found this approach to work well (using software):
Find the midplane of 3D space between direction vectors [ 1 0 0 ] and newXaxis.
Find the midplane of 3D space between direction vectors [ 0 1 0 ] and newYaxis.
The intersection of the two midplanes (a 3D line) gives the Euler axis direction.
Construct a plane thru the origin orthogonal to the Euler axis. Project any of the "before-after" pairs onto direction vectors on this plane, e.g. [ 1 0 0 ] and newXaxis. The cos(Euler angle) is the dot product of these 2 projected direction vectors. The sin(Euler angle) is the magnitude of their cross product. You can use an IEEE "atan2(sin, cos)" function to obtain an accurate measure of the Euler angle.
If algorithmizing this conversion, be sure to pick off the edge cases before plowing forward into the midplane calcs. These cases are the identity rotator, and the 3 single axis rotators (e.g., rotate about x-axis, about negative x-axis, y-axis, negative y-axis, etc.)