Last weekend I searched for the reason why skew and quaternion results were different in my methods. The reason was in the end simple, as always in retrospect... The skew symmetric matrix must be normalized for the published formulae to work correctly (the w=1).
Now there is a nice match between quaternions and skew and I can proceed.
Let me describe how my development works, if someone is interested:
I try to find test data in the internet, but in most cases without success. So I have to build them my own, which is very time consuming. Life would be much easier if for formulae an example would be given.
I create methods, rotationToQuaternions for example, and the oppsosite, quaternionsToRotations (I call it forward and inverse methods)
I create random data/matrices/angles/properties million of times and verify that applying the forward and inverse methods result in the starting data. Rounding errors must be accepted, but in most cases if the errors are too high, it is due to a bug in the methods. For float, I accept errors in the region of 1e-7 to 1e-8 for single calculations and worst 1e-4 if the calculation is very complex. A very common bug is confusing radians with sin or cos values, because they are similar for low angles.
a different validation method is to compare my implementation with a library, I e. g. compared the results of the generalized inverse with Eigen
another test method is unit testing in the spirit of JUnit. I use it when refactoring methods as well
Unfortunately, there are books with errors in the formulae (I have one book with an error on every third page), so I have to verify the formulae. With examples this could be deteced.
An enlithing article to understand the difference between ZYX transformations and infinitesimal angle changes (DH/rotation matrix vs quaternion/skew) is "Derivative of Rotation Matrix... by Fumio Hamano. I can recommend reading it. Both aspects build the core of the robot kinematics.
I have to clean up the code again and verify that the kinematics works.