Yes sure. see what the constraint nodes do, they take the parent matrix and the individual channels in for the job, this ensures the plugs are clean.
Secondly, I would have to accumulate all the rotations along DAG manually - which seems like a really bad idea.
Well that may or may not be possible it really depends what you need. Thing is rotation as such does not exist, the rotation is off course encoded in the matrix as orientation which for all intents and purposes EXCEPT interpolation is the same thing. This however exposes you to a big logic hole, see rotation does not exist really at all. Which is why you cant really compute anything on top of euler rotations, save for situations where only one channel changes. Reason being there's not 3 axes of rotation there is just one rotation ever present. Now matrices and quaternions work right but they dont measure rotation, they measure orientation. Also note even mr Euler didn't CLAIM anything other than describing orientation.
This problem can not be generally solved, i suspect because it cant but can not prove that off course. Tough you can find a lot of special solutions that work in limited scope. So either you need to find your own set of mathematical equations* or orientation is enough for you in which case you can extract the rotation form the matrix.
which seems like a really bad idea.
Maya does in this all the time in form of matrices.
- and then your on your own, uncharted scientific area. Publish your stuff on journals it futhers science.