//Author: Daniel Hughes

function Clifford(d)
FF := CyclotomicField(LeastCommonMultiple(8*3,2*d));
ZZ := IntegerRing();
// The factor of 8 ensures sqrt(d) is in the field
// The factor of 3 ensures the third roots of unity are
PR<x>:=PolynomialRing(FF);
rt:=Roots(x^2-d)[1][1];
// make sure rt is indeed the positive square root of d
if Real(rt) lt 0 then rt := -rt; end if;

w:=RootOfUnity(d); mu:=RootOfUnity(2*d);
zeta:=RootOfUnity(24);

I := Matrix(FF,d,d,[<j+1,j+1,1>: j in [0..d-1]]);
A := Matrix(FF,d,d,[<j+1,j+1,w^j>: j in [0..d-1]]);
S := Matrix(FF,d,d,[<((j+1) mod d)+1,j+1,1>: j in [0..d-1]]);
F := Matrix(FF,d,d,[<j+1,k+1,w^(j*k)>: j,k in [0..d-1]])/rt;
R := Matrix(FF,d,d,[<j+1,j+1,mu^(j*(j+d))>: j in [0..d-1]]);
M := R*F;  Z := zeta^(d-1)*M;

G:=MatrixGroup< d,FF|F,R>;
// This is the Clifford group
G:=MatrixGroup< d,FF|S,A,F,R>;
return G;
end function;

function tClifford(c,i)
P:=Clifford(c);
C:=Centre(P);
T:=RightTransversal(P,C);
a:=Basis(Eigenspace(Identity(P),1))[i];
v:=[a^g:g in T];
V:=Parent(v[1]); subs:={@ @}; vecs:=[ ];
for j:=1 to #v do
S:=sub<V|v[j]>;
if S notin subs then
Include(~subs,S);Append(~vecs,v[j]);
end if;
end for;
W:=Transpose(Matrix(vecs));
A:=CanonicalGramian(W);
n:=Ncols(W);
d:=Nrows(W);

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(A) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else
c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);
end if;
n:=Ncols(A); sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (A[i,j]*ComplexConjugate(A[i,j]))^(t);
end for;
end for;

D:=Diagonal(A); diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + D[i]^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;

end while;
return t-1,n;
end function;

function Matip(A,B)
ip:=Trace(A*HermTranspose(B));
return ip;
end function;


function tCliffordMat(c);
P:=Clifford(c);
C:=Centre(P);
T:=RightTransversal(P,C);
d:=c^2;

dif:=0; t:=0;
while dif eq 0 do
t:=t+1;
if BaseRing(T[1]) eq Rationals() then
c:=&*[(1+i)/(d+i): i in [0..2*(t-1) by 2]];
else
c:=Factorial(d-1)*Factorial(t)/Factorial(d+t-1);
end if;

n:=#T; sum:=0;
for i:=1 to n do
for j:=1 to n do
sum:=sum + (Matip(T[i],T[j])*ComplexConjugate(Matip(T[i],T[j])))^(t);
end for;
end for;

diagsum:=0;
for i:=1 to n do
diagsum:=diagsum + Matip(T[i],T[i])^t;
end for;
diagsum:=diagsum^2;

dif:=sum-c*diagsum;
end while;
return t-1,n;
end function;