Interpolating 3D matrices, quaternions?

Had anyone seen some free quaternion based code doing thing
similar to player v10 Matrix3D interpolate?

In my tests, based on pv3d Matrix3D.rotationMatrix(x, y, z, angle) or
similar, there is significant divergency of basis vectors. Simple
orthonormalization a-la Gram/Schmidt prevents divergency, but adds
"preferred" direction into cocktail, and other suggested methods of
orthonormalization
(http://sci.tech-archive.net/Archive/.../msg00114.html)
are getting too complex for this should-have-been-simple thing. So, I
thought, maybe quaternion implementation will not suffer this problem?
This is not engine-specific, so please post whatever comes close.