printf(" LinearAlgebra, quat(a,b,c,d), qconj(x),qexp(u,theta), qmul(p,q), qMul(q1, q2,...,qn).\n");


#Note that the following matrix representation of a quaternion is not unique. See wikipedia quaternion matrix representation.

with(LinearAlgebra):

quat:=proc(a,b,c,d)
Matrix([
[a,b,c,d],
[-b,a,-d,c],
[-c,d,a,-b],
[-d,-c, b,a]
]);
end proc:

qconj:=proc(x) 
[x[1],-x[2],-x[3],-x[4]];
end:

qexp:=proc(u,theta)
[cos(theta), sin(theta)*u[1],sin(theta)*u[2],sin(theta)*u[3]];
end proc:

qmul:=proc(A,X) 
local q;
description "multiplies two quaternions";
q:=quat(op(A)).quat(op(X));
[q[1,1],q[1,2],q[1,3],q[1,4]];
map(simplify,%);
end proc:

qMul:=proc()
local q,i;
description "multiplies finite number of quaternions";
q:=[1,0,0,0];
for i from 1 to _npassed do
q:=qmul(q,_passed[i]);
od:
q;
end proc: