###############################################################
#
#      Multiplication rules for quaternions, 
#          Cayley numbers, and octonians 
#                 for MAPLE V5
#
# References: chapter 22 of Fulton and Harris, 
# Representation Theory, Springer-Verlag, 1991.
# page 45 of R. Schafer, 
# An introduction to non-associative algebras, Dover, 1994.
# section 6.3 of I. Kantor and A. Solodovnikov, 
# Hypercomplex numbers, Springer-Verlag, 1989.
#
# wdj, cayley.mpl   
###############################################################
#
# A Cayley number will be written as a 2x2 matrix
# with real numbers on the diagonal and 3-vectors
# on the skew-diagonal in the form
# x:=[x1,x2,x3];y:=[y1,y2,y3];
# v:=array(1..2,1..2,[[a,x],[y,b]]);
# where a,b,x1,...,y3 are real

multiply_cayley:=proc(A::array,B::array)
#assumes arrays are Cayley numbers
 local a1,a2,b1,b2,x1,x2,y1,y2,C,
 x11,x12,x13,y11,y12,y13,x21,x22,x23,y21,y22,y23;
 a1:=A[1,1];
 b1:=A[2,2];
 x1:=A[2,1];
 y1:=A[1,2];
 a2:=B[1,1];
 b2:=B[2,2];
 x2:=B[2,1];
 y2:=B[1,2];
 C:=array(1..2,1..2,[[a1*a2+dotprod(x1,y2,`orthogonal`),expand(a1*x2+b2*x1)],[expand(a2*x1+b1*x2),b1*b2+dotprod(x2,y1,`orthogonal`)]]);
 RETURN(evalm(C));
end:

cayley_basis:=proc()
 local e1,e2,e3,e4,e5,e6,e7,e8;
 e1:=array(1..2,1..2,[[0,[1,0,0]],[[0,0,0],0]]);
 e2:=array(1..2,1..2,[[0,[0,1,0]],[[0,0,0],0]]);
 e3:=array(1..2,1..2,[[0,[0,0,1]],[[0,0,0],0]]);
 e4:=array(1..2,1..2,[[0,[0,0,0]],[[1,0,0],0]]);
 e5:=array(1..2,1..2,[[0,[0,0,0]],[[0,1,0],0]]);
 e6:=array(1..2,1..2,[[0,[0,0,0]],[[0,0,1],0]]);
 e7:=array(1..2,1..2,[[1,[0,0,0]],[[0,0,0],0]]);
 e8:=array(1..2,1..2,[[0,[0,0,0]],[[0,0,0],1]]);
 RETURN([evalm(e1),evalm(e2),evalm(e3),evalm(e4),evalm(e5),evalm(e6),evalm(e7),evalm(e8)]);
end:

cayley_multiplication_table:=
array(1..9,1..9,[[C,seq(e[i],i=1..8)],
 [e[1],0,0,0,e[8],0,0,0,0],
 [e[2],0,0,0,0,e[8],0,0,0],
 [e[3],0,0,0,0,0,e[8],0,0],
 [e[4],e[7],0,0,0,0,0,e[4],e[1]],
 [e[5],0,e[7],0,0,0,0,e[5],e[1]],
 [e[6],0,0,e[7],0,0,0,e[6],e[2]],
 [e[7],0,0,0,e[1],e[2],e[3],e[7],0],
 [e[8],0,0,0,e[4],e[5],e[6],0,e[8]]]):

conjugate_cayley:=proc(w::array)
#assumes array is a Cayley number
 local w0,x,y,a,b;
 a:=w[1,1]; b:=w[2,2];
 x:=w[1,2]; y:=w[2,1];
 w0:=array(1..2,1..2,[[a,y],[x,b]]);
 RETURN(evalm(w0));
end:

norm_cayley:=proc(w::array)
#assumes array is a Cayley number
 local w0,w1;
 w0:=conjugate_cayley(w);
 w1:=multiply_cayley(w,w0);
 RETURN(trace(w1));
end:

trace_cayley:=proc(w::array)
#assumes array is a Cayley number
 local w0,w1;
 w0:=conjugate_cayley(w);
 w1:=(w+w0)/2;
 RETURN(evalm(w1));
end:

form_cayley:=proc(w1::array,w2::array)
#assumes arrays are Cayley numbers
 local a,b,form;
 a:=multiply_cayley(w1,conjugate_cayley(w2));
 b:=multiply_cayley(w2,conjugate_cayley(w1));
 form:=(a+b)/2;
 RETURN(evalm(form));
end:

#
#the bilinear form below is the matrix trace
#of the form above
#
#G_2 is the group of *algebra* automorphisms
#of the cayley numbers preserving the following form.
#
#SO(8) is the group of *vector space* automorphisms
#of the cayley numbers preserving this form.
#
#the cayley_basis command above returns an orthonormal
#basis with respect to this form
#

bilinearform_cayley:=proc(w1::array,w2::array)
#assumes arrays are Cayley numbers
 local x1,y1,a1,b1,x2,y2,a2,b2,form;
 a1:=w1[1,1]; b1:=w1[2,2];
 x1:=w1[1,2]; y1:=w1[2,1];
 a2:=w2[1,1]; b2:=w2[2,2];
 x2:=w2[1,2]; y2:=w2[2,1];
 form:=a1*a2+b1*b2+dotprod(x1,x2,`orthogonal`)+dotprod(y1,y2,`orthogonal`);
 RETURN(form);
end:

conjugate_complex:=proc(w)
#assumes w is a complex number x+I*y
 local w0,x,y,r,s,t,u;
 s:=coeff(w,I); r:=w-I*s;
 w0:=r-s*I;
 RETURN(w0);
end:

# A quaternion will be written as 
# w:=[r+I*s,t+I*u];
# regarded as a+j*b, where a and b are complex and 
# j is the usual unit quaternion j.

multiply_quaternion:=proc(w1::list,w2::list)
#assumes lists w1, w2 are quaternion numbers
 local r1,s1,t1,u1,r2,s2,t2,u2,w,x1,x2,y1,y2;
 x1:=w1[1];y1:=w1[2];x2:=w2[1]; y2:=w2[2];
 w:=[evalc(expand(x1*x2-y1*conjugate_complex(y2))),evalc(expand(conjugate_complex(x2)*y1+y2*x1))];
 RETURN(w);
end:

# The unit quaternions are
q1:=[1,0]: qi:=[I,0]: qj:=[0,1]: qk:=[0,I]:

quaternion_basis:=proc()
 local e1,e2,e3,e4,e5,e6,e7,e8;
 e1:=[1,0];
 e2:=[I,0];
 e3:=[0,1];
 e4:=[0,I];
 RETURN([e1,e2,e3,e4]);
end:

conjugate_quaternion:=proc(w::list)
#assumes list is a quaternion number
 local w0,x,y,r,s,t,u;
 x:=w[1];s:=coeff(x,I); r:=x-I*s;
 y:=w[2];u:=coeff(y,I); t:=y-I*u;
 w0:=[r-s*I,-t-u*I];
 RETURN(w0);
end:

norm_quaternion:=proc(w::list)
#assumes list is a quaternion number
 local w0,w1;
 w0:=conjugate_quaternion(w);
 w1:=multiply_quaternion(w,w0);
 RETURN(op(1,w1));
end:

trace_quaternion:=proc(w::list)
#assumes list is a quaternion number
 local w0,w1;
 w0:=conjugate_quaternion(w);
 w1:=(w+w0)/2;
 RETURN(op(1,w1));
end:

# An octonian number will be represented in the form
# o:=[[r1+I*s1,t1+I*u1],[r2+I*s2,t2+I*u2]];
# where r1,s1,r2,s2,t1,t2,u1,u2 are real.

multiply_octonian:=proc(oo1::list,oo2::list)
#assumes lists are octonian numbers
#uses muliplication rule in FH and KS
#(the rule in Schafer appears differently)
 local a,b,c,d,w;
 a:=oo1[1];b:=oo1[2]; c:=oo2[1];d:=oo2[2];
 w:=[evalc(expand(multiply_quaternion(a,c)-multiply_quaternion(conjugate_quaternion(d),b))),evalc(expand(multiply_quaternion(d,a)+multiply_quaternion(b,conjugate_quaternion(c))))];
 RETURN(w);
end:

octonian_basis:=proc()
 local e1,e2,e3,e4,e5,e6,e7,e8;
 e1:=[[1,0],[0,0]];
 e2:=[[I,0],[0,0]];
 e3:=[[0,1],[0,0]];
 e4:=[[0,I],[0,0]];
 e5:=[[0,0],[1,0]];
 e6:=[[0,0],[I,0]];
 e7:=[[0,0],[0,1]];
 e8:=[[0,0],[0,I]];
 RETURN([e1,e2,e3,e4,e5,e6,e7,e8]);
end:

conjugate_octonian:=proc(oo::list)
#assumes list is octonian number
 local a,b,w;
 a:=oo[1];b:=oo[2]; w:=[conjugate_quaternion(a),-b];
 RETURN(w);
end:

norm_octonian:=proc(w::list)
#assumes list is a octonian number
 local w0,w1;
 w0:=conjugate_octonian(w);
 w1:=multiply_octonian(w,w0);
 RETURN(op(1,op(1,w1)));
end:

trace_octonian:=proc(w::list)
#assumes list is a octonian number
 local w0,w1;
 w0:=conjugate_octonian(w);
 w1:=(w+w0)/2;
 RETURN(op(1,op(1,w1)));
end: