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

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

kk:=Rationals;

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

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

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

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

phiZm^3-3*phiZm-3 = Zero(mm);    
phi2Zm^3-3*phi2Zm-3 = Zero(mm);    

phiZm =  (2/5*Sqrt(-15))*One(mm)
        +(-1/2+3/10*Sqrt(-15))*Zm
        +(-1/5*Sqrt(-15))*Zm^2;
phi2Zm = (-2/5*Sqrt(-15))*One(mm)
        +(-1/2-3/10*Sqrt(-15))*Zm
        +(1/5*Sqrt(-15))*Zm^2;

phi2Zm = (2/5*Sqrt(-15))*One(mm)
        +(-1/2+3/10*Sqrt(-15))*phiZm
        +(-1/5*Sqrt(-15))*phiZm^2;
Zm =     (2/5*Sqrt(-15))*One(mm)
        +(-1/2+3/10*Sqrt(-15))*phi2Zm
        +(-1/5*Sqrt(-15))*phi2Zm^2;
phiZm =  (-2/5*Sqrt(-15))*One(mm)
        +(-1/2-3/10*Sqrt(-15))*phi2Zm
        +(1/5*Sqrt(-15))*phi2Zm^2;
Zm =     (-2/5*Sqrt(-15))*One(mm)
        +(-1/2-3/10*Sqrt(-15))*phiZm
        +(1/5*Sqrt(-15))*phiZm^2;

# 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:=  (2/5*Sqrt(-15))*I3
         +((-1/2+3/10*Sqrt(-15)))*ZDM
         +((-1/5*Sqrt(-15)))*ZDM^2;
phi2ZDM:= ((-2/5*Sqrt(-15)))*I3
         +(-1/2-3/10*Sqrt(-15))*ZDM
         +(1/5*Sqrt(-15))*ZDM^2;

ZDM^3-3*ZDM-3*One(ZDM) = Zero(ZDM);
phiZDM^3-3*phiZDM-3*One(ZDM) = Zero(ZDM);
phi2ZDM^3-3*phi2ZDM-3*One(ZDM) = Zero(ZDM) ;

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

sigmaDM^3 = 2*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:= [[2,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(-15)*One(mm)) = -I3*Sqrt(-15)*One(mm);

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

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

VecToDM:=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(-15)^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);