A different algorithm that offers the same functionality has already
been implemented earlier by David Kohel, see the Quaternion Algebras
package on:

http://www.maths.usyd.edu.au/u/kohel/alg/index.html