q1={4,{2,0,-1}}
q2={1,{-2,3,1}}

Qprod[q1_,q2_]:={q1[[1]]*q2[[1]]-q1[[2]].q2[[2]], Cross[q1[[2]],q2[[2]]]+q1[[1]]*q2[[2]]+q2[[1]]*q1[[2]]}

Qprod[q1,q2]

Qnorm[q_]:=Sqrt[q[[1]]^2+q[[2]].q[[2]]]

Qnorm[q1]

Qcon[q_]:={q[[1]],-q[[2]]}

Qcon[q1]

Qrot[q_,v_]:=Qprod[Qprod[q,v],Qcon[q]]

v={0,{1,2,3}}

Qrot[q1,v]

q={s,{x,y,z}}
q=q/Qnorm[q]

Qprod[q1,q2]

QtoMatrix[q_]:={{s,-x,z,-y},{x,s,y,z},{-z,-y,s,x},{y,-z,-x,s}}/.{s->q[[1]], x->q[[2]][[1]], y->q[[2]][[2]], z->q[[2]][[3]]}

QtoMatrix[q1]//MatrixForm

QtoMatrix[q2].QtoMatrix[q1]

A=Transpose[{Simplify[Qprod[Qprod[q,{0,{1,0,0}}], Con[q]]][[2]], Simplify[Qprod[Qprod[q,{0,{0,1,0}}], Con[q]]][[2]], Simplify[Qprod[Qprod[q,{0,{0,0,1}}], Con[q]]][[2]]}]

MatrixForm[A]