#####
#
# This is GAP code
# 
# The file is ~steger/donald/SU21/a1p5/a1p5-proj.gap
#
# See file ~steger/donald/SU21/a1p5/D.txt for the definition of the
# division algebra $\D$ and for all the notation.
#
#####

#####
#
# $k=\QQ
#

kk:=Rationals;

#####
#
# $\ell=\QQ[\sqrt{-2}]$
#
ll:=Field(Sqrt(-1));
llBasis:=Basis(ll,[1,Sqrt(-1)]);

#####
#
# $m=\ell[Z]/(Z^3-3*Z^2-2)
#
xl:=Indeterminate(ll,"xl");
f:=xl^3-3*xl^2-2*One(ll); # Discriminant = -324 = -18^2
mm:=FieldExtension(ll,f);
Zm:=RootOfDefiningPolynomial(mm);

mmBasis:=Basis(mm,[One(mm),Zm,Zm^2]);

fFcn:=x->x^3-3*x^2-2*One(x);
fFcn(Zm) = Zero(mm);

mmCoef:=x->Concatenation(
  List(Coefficients(mmBasis,x),y->Coefficients(llBasis,y)));

mmLinComb:=vec->LinearCombination(mmBasis,[
  LinearCombination(llBasis,[vec[1],vec[2]]),
  LinearCombination(llBasis,[vec[3],vec[4]]),
  LinearCombination(llBasis,[vec[5],vec[6]])
]);

#####
#
# $\phi:m\to m$
#

sum:=3-Zm; # Sum of the other two roots
diff:=Sqrt(Discriminant(f))/(Zm^2-sum*Zm+2/Zm);
  # Difference of the other two roots
phiZm:=(sum-diff)/2;
phi2Zm:=(sum+diff)/2;

phiOnZ:=x->
  (3/2+1/2*Sqrt(-1))*One(x)
        +(-1/2-2*Sqrt(-1))*x
        +(1/2*Sqrt(-1))*x^2;

phi2OnZ:=x->
  (3/2-1/2*Sqrt(-1))*One(x)
        +(-1/2+2*Sqrt(-1))*x
        +(-1/2*Sqrt(-1))*x^2;

IsZero(fFcn(Zm));
IsZero(fFcn(phiZm));
IsZero(fFcn(phi2Zm));

phiOnZ(Zm) = phiZm;
phiOnZ(phiZm) = phi2Zm;
phiOnZ(phi2Zm) = Zm;

phi2OnZ(Zm) = phi2Zm;
phi2OnZ(phi2Zm) = phiZm;
phi2OnZ(phiZm) = Zm;

# Convert an element which belongs of~$m$ which actually belongs
# to~$\ell$ to the corresponding element of~$\ell$
mTol:=function(x)
  local c;
  c:=Coefficients(mmBasis,x);
  if c[2]<>0 or c[3]<>0 then
    return fail;
  else
    return c[1];
  fi;
end;

#####
#
# Conjugation in~$m$ (taking $\bar Z=Z$)
#
mmConj:=function(x)
  local c;
  c:=ExtRepOfObj(x);
  return ComplexConjugate(c[1])*One(mm)
        +ComplexConjugate(c[2])*Zm
        +ComplexConjugate(c[3])*Zm^2;
end;

mmConj(Zm) = Zm;
mmConj(phiZm) = phi2Zm;
mmConj(phi2Zm) = phiZm;

#####
#
# Utility functions for matrices
#

I3:=IdentityMat(3); # 3 by 3 identity matrix

IsPOne:=mtx->mtx=mtx[1][1]*One(mtx);

fOnM:=function(f,mtx)
  return List(mtx,row->List(row,x->f(x)));
end;

# Function to calculate the adjoint of a matrix over~$m$
mmAdj:=mtx->TransposedMat(List(mtx,row->List(row,x->mmConj(x))));

######
#
# $Z$, $\phi(Z)$, and $\phi^2(Z)$ in the matrix representation
# of~$\D$ over~$m$
#
ZDM:=[[Zm,0,0],[0,phiZm,0],[0,0,phi2Zm]]*One(mm);
phiZDM:=phiOnZ(ZDM);
phi2ZDM:=phi2OnZ(ZDM);

IsZero(fFcn(ZDM));
IsZero(fFcn(phiZDM));
IsZero(fFcn(phi2ZDM));

#####
#
# $\sigma$ in the matrix representation of~$\D$ over~$m$
#
sigmaDM:=[[0,1,0],[0,0,1],[5,0,0]]*One(mm);

sigmaDM^3 = 5*One(ZDM);
sigmaDM*ZDM = phiZDM*sigmaDM;
sigmaDM*phiZDM = phi2ZDM*sigmaDM;
sigmaDM*phi2ZDM = ZDM*sigmaDM;

#####
#
# In the matrix representation of~$\D$ over~$m$, the matrix~$F$ so
# that
#
#   \iota(A) = F^{-1} A^* F
#
FDM:=[[5,0,0],[0,0,1],[0,1,0]]*One(mm);

iotaDM:=mtx->FDM^-1*mmAdj(mtx)*FDM;

mmAdj(FDM) = FDM;
iotaDM(sigmaDM) = sigmaDM;
iotaDM(ZDM) = ZDM;
iotaDM(phiZDM) = phi2ZDM;
iotaDM(phi2ZDM) = phiZDM;
iotaDM(I3*Sqrt(-1)*One(mm)) = -I3*Sqrt(-1)*One(mm);

#####
#
# Set up $\D$
#
#####

DMtoDMVec:=mtx->Concatenation(List(mtx,
  row->Concatenation(List(row,x->mmCoef(x)))));

DMVecToDM:=vec->[
  [mmLinComb(vec{[1..6]}),
     mmLinComb(vec{[7..12]}),
       mmLinComb(vec{[13..18]})], 
  [mmLinComb(vec{[19..24]}),
     mmLinComb(vec{[25..30]}),
       mmLinComb(vec{[31..36]})], 
  [mmLinComb(vec{[37..42]}),
     mmLinComb(vec{[43..48]}),
        mmLinComb(vec{[49..54]})] 
];

# 18-element basis of $\D$ over $\QQ$
basisDM:=ListX([0..1],[0..2],[0..2],
  function(j,k,l)
    return Sqrt(-1)^j*ZDM^k*sigmaDM^l;
  end
);;

basisDMVec0:=List(basisDM,DMtoDMVec);;
DMVec:=VectorSpace(Rationals,basisDMVec0);
basisDMVec:=Basis(DMVec,basisDMVec0);;

#####
#
# Find the matrix for $\iota$ acting on $\D$
#
#####

iotaPDM:=TransposedMat(
  List(basisDM,bmtx->
    Coefficients(basisDMVec,DMtoDMVec(iotaDM(bmtx)))
));
	
IsOne(iotaPDM^2);

######
#
# Contruct the the self-adjoint, traceless part of~$\D$
#
#####

TracesDM:=List(basisDM,mtx->Trace(ll,kk,mTol(TraceMat(mtx))));

DMVecSelfAdjBasis0:=NullspaceMat(TransposedMat(
  Concatenation(One(iotaPDM)-iotaPDM,[TracesDM])));;
DMVecSelfAdj:=VectorSpace(Rationals,DMVecSelfAdjBasis0);
DMVecSelfAdjBasis:=Basis(DMVecSelfAdj);