I think I found a bug in Quaternion:
Before calculating the rotation matrix from an quaternion, you need to normalize it first. It looks like in both implementations of Quaternion.toRotationMatrix this is attempted by dividing (a little bit indirectly) by Quaternion.norm(). This seems to be not correct, as norm() is defined as w