RecToList:=function(A); //Given the record that Subgroups or MaximalSubgroups returns, this returns a list of the subgroups, rather than all the stuff the record has.
L:=[];
for i:=1 to #A do
L:=Append(L,A[i]`subgroup);
end for;
return L;
end function;

ReadInFile:=function(s); //This lets me read in files easily.
k:=Read(s);
return eval k;
end function;

conj:=function(h,g); //This just conjugates h by g and is shorthand for me.
return g*h*g^-1;
end function;

ListOfIdentities:=function(n) //Create a list with n 2x2 identity matrices in it. This will be useful in later functions.
I:=IdentityMatrix(Integers(),2);
N:=[];
for i:=1 to n do
N:=Append(N,I);
end for;
return N;
end function;

TensorProductList:=function(N); //This does the tensor product of the matrices in the list N.
T:=N[1];
for k:=2 to #N do
T:=TensorProduct(T,N[k]);
end for;
return T;
end function;

sigma:=function(i,j,m); //This gives an m-fold tensor product with little sigma_i in its jth slot and a 2x2 identity in every other slot.
if i eq 1 then
s:=Matrix([[0,1],[1,0]]);
else
s:=Matrix([[1,0],[0,-1]]);
end if;
N:=ListOfIdentities(m);
N[j]:=s;
return ChangeRing(TensorProductList(N),Rationals());
end function;

E:=function(m); //This gives the fundamental representation of the extraspecial 2-group of degree m by using the sigma function to compute each of the 2m generators.
Gens:=[];
for a:=1 to m do
for b:=1 to m do
Gens:=Append(Gens,sigma(1,a,m));
Gens:=Append(Gens,sigma(2,b,m));
end for;
end for;
return sub<GeneralLinearGroup(2^m,Rationals())|Gens>;
end function;

g1:=function(m); //This computes the first new generator of BRW(m). It's a block diagonal matrix with an identity matrix, a zero, a zero, and (I x I x ... x sigma_1).
if m eq 1 then
return IdentityMatrix(Integers(),2);
end if;
I:=IdentityMatrix(Integers(),2^m div 2);
Z:=ZeroMatrix(Integers(),2^m div 2);
N:=ListOfIdentities(m-1);
N[m-1]:=sigma(1,1,1);
T:=TensorProductList(N);
return BlockMatrix(2,2,[I,Z,Z,T]);
end function;

g2:=function(m); // This computes the generator corresponding to the second new generator
p:=[1..2^m]; //Initialize permutation to use to compute the generator
for i:=1 to 2^(m-1) do //Creates the permutation on 2^m elements corresponding to the second generator of GL(m,2)'s action on R^(2^m), it is [1,3,5,...,2^(m-1),2,4,...,2^m]
p[i]:=2*i-1;
p[i+2^(m-1)]:=2*i;
end for;
return Transpose(PermutationMatrix(Rationals(),p)); //Returns the corresponding matrix representation of this permutation
end function;

D:=function(m); //This gives a 2^m identity matrix, except each 4th entry is a -1.
d:=IdentityMatrix(Integers(),2^m);
for i:=1 to Floor(2^m/4) do
d[4*i][4*i]:=-1;
end for;
return d;
end function;

H:=function(m) //This creates the Hadamard matrix.
Q<w>:=QuadraticField(2);
h:=1/w*Matrix([[1,1],[1,-1]]);
for i:=1 to m-1 do //Creates the Kronecker product of h with I m-1 times
h:=KroneckerProduct(IdentityMatrix(Q,2),h);
end for;
return h;
end function;

BRW:=function(m); //This creates the BRW group. It is generated by the gens of E(m); and g1, g2, and D, plus each of those conjugated by the Hadamard matrix.
return sub<GeneralLinearGroup(2^m,Rationals())|Generators(E(m)),g1(m),g2(m),D(m),conj(g1(m),H(m)),conj(g2(m),H(m)),conj(D(m),H(m))>;
end function;

CG:=function(m); //This creates the Clifford group. It is generated by the gens of BRW and the Hadamard matrix.
return sub<GeneralLinearGroup(2^m,QuadraticField(2))|Generators(BRW(m)),H(m)>;
end function;

vacuumsum:=function(M); //This computes the sum part of the formula for the vacuum anomaly of M.
M:=Matrix(M);
n:=Order(M);
D:=0;
if n ne 1 then
D:=1;
end if;
S:=0;
if D ne 0 then
if n ne 2 then
C:=CyclotomicPolynomial(n);
F<p>:=SplittingField([C]);
M:=ChangeRing(M,F);
else
p:=-1;
end if;
S:=0;
for k:=1 to n-1 do
S:=S+k*(n-k)*Dimension(Eigenspace(M,p^k));
end for;
end if;
return S;
end function;

vacuum:=function(M); //This computes the vacuum anomaly of M.
M:=Matrix(M);
S:=vacuumsum(M);
n:=Order(M);
return 1/(4*n^2) * vacuumsum(M);
end function;

type:=function(M); //This computes the type of M.
M:=Matrix(M);
V:=vacuum(M);
n:=Order(M);
k:=n^2*V;
if IsCoercible(Integers(),k) then
l:=(Integers() ! k) mod n;
else
l:=-1;
end if;
return l;
end function;

chi1:=function(g); //This is shorthand for the character of the fundamental representation.
return Trace(Matrix(g));
end function;

chi2:=function(g); //This computes the sum part of the inner product of the trivial character with the symmetric character.
M:=Matrix(g);
L:=M^2;
a:=Trace(M);
b:=Trace(M^2);
c:=Trace(M)^2;
return a+(b+c)/2;
end function;

dimchi1:=function(H); //This computes the inner product of the fundamental character with the trivial character.
C:=ConjugacyClasses(H);
z:=0;
for i:=1 to # C do
a:=chi1(C[i][3]);
z:=z+a*C[i][2];
end for;
return z/Order(H);
end function;

dimchi2:=function(H); //This computes the inner product of the trivial character with the symmetric character.
C:=ConjugacyClasses(H);
z:=0;
for k:=1 to # C do
a:=chi2(C[k][3]);
z:=z+a*C[k][2];
end for;
return z/Order(H);
end function;

CheckConjugacyClasses:=function(H); //This checks to see if H, a group, has any elements with vacuum anomaly <= number. (I manually adjusted number as I needed it.)
number:=1;
for h in Generators(H) do //Check the relatively small number of generators first before generating all the conjugacy classes!
if vacuum(h) le number then
return false;
end if;
end for;
C:=ConjugacyClasses(H);
for i:=2 to # C do
if vacuum(C[i][3]) le number then
return false;
end if;
end for;
return true;
end function;

CheckDimensions:=function(H); //This checks to see if the inner products of the fundamental with the trivial and the symmetric with the trivial are 0 and 1, respectively.
if dimchi1(H) ne 0 then
return false;
end if;
if dimchi2(H) ne 1 then
return false;
end if;
return true;
end function;

MatrixToCycle:=function(M); //Elements of E(m) are signed permutation matrices and can hence be represented as permutations of 1 up to 2^(m+1). The first half represents the columns, and the second represents the minu of the columns.
n:=NumberOfRows(M);
Z:=ZeroMatrix(Integers(),n,1);
N:=[];
for i:=1 to NumberOfRows(Z) do
k:=Z;
k[i,1]:=1;
N:=Append(N,k);
end for;
for i:=1 to NumberOfRows(Z) do
k:=Z;
k[i,1]:=-1;
N:=Append(N,k);
end for;
A:=[];
M:=ChangeRing(M,Integers());
K:=[];
for i:=1 to #N do
k:=M*N[i];
for j:=1 to #N do
if k eq N[j] then
K:=Append(K,j);
end if;
end for;
end for;
return SymmetricGroup(n*2) ! K;
end function;

CycleToMatrix:=function(A); //This goes the other way if we ever need to.
n:=Degree(Parent(A)) div 2;
Z:=ZeroMatrix(Integers(),n,1);
N:=[];
for i:=1 to NumberOfRows(Z) do
k:=Z;
k[i,1]:=1;
N:=Append(N,k);
end for;
for i:=1 to NumberOfRows(Z) do
k:=Z;
k[i,1]:=-1;
N:=Append(N,k);
end for;
M:=[];
for i:=1 to n do
M:=Append(M,N[i ^ A]);
end for;
return Transpose(Matrix(M));
end function;

PermE:=function(m); //Returns the permutation version of E(m). Calculations with this run about 2.5 times faster than with matrices.
Gens:=[];
for x in Generators(E(m)) do
Gens:=Append(Gens,MatrixToCycle(x));
end for;
return sub<SymmetricGroup(2^m*2)|Gens>;
end function;

MyE:=function(m); //This is for my convenience since I relabeled the function.
return PermE(m);
end function;

FindTrace:=function(A); //This finds the trace of a permutation.
z:=0;
n:=Degree(Parent(A)) div 2;
for i:=1 to n do
if i ^ A eq i then
z:=z+1;
end if;
if i ^ A eq (i + n) then
z:=z-1;
end if;
end for;
return z;
end function;

chi2perm:=function(g); //Computes the inner product of the symmetric with the trivial for the permutation realization.
a:=FindTrace(g);
b:=FindTrace(g^2);
c:=FindTrace(g)^2;
return a+(b+c)/2;
end function;

dimchi1perm:=function(H); //Computes the inner product of the fundamental with the trivial for the permutation realization.
C:=ConjugacyClasses(H);
z:=0;
for i:=1 to # C do
a:=FindTrace(C[i][3]);
z:=z+a*C[i][2];
end for;
return z/Order(H);
end function;

dimchi2perm:=function(H); //Computes the inner product of the symmetric with the trivial for the permutation realization.
C:=ConjugacyClasses(H);
z:=0;
for k:=1 to # C do
a:=chi2perm(C[k][3]);
z:=z+a*C[k][2];
end for;
return z/Order(H);
end function;

TranslateSubgroupOfCycles:=function(H); //This takes a subgroup in the permutation realization to the corresponding subgroup in the matrix realization.
m:=Integers() ! Log(2,Degree(Parent(Random(H))) div 2);
L:=[];
for c in Generators(H) do
L:=Append(L,CycleToMatrix(c));
end for;
return sub<E(m)|L>;
end function;

balanced:=function(m); //This returns the generator matrix for the Balanced Barnes-Wall lattice (as rows).
Q<w>:=QuadraticField(2);
M:=Matrix([[w,0],[1,1]]);
T:=M;
for i:=1 to m-1 do//
T:=TensorProduct(T,M);
end for;
return T;
end function;

Nebegen:=function(m); //This returns the generator matrix for Nebe et al's construction of the Barnes-Wall lattice.
Q<w>:=QuadraticField(2);
T:=balanced(m);
B:=T;
for i:=1 to 2^m do
for j:=1 to 2^m do
if Eltseq(T[i,j])[2] ne 0 then //Eltseq[2] tells me the sqrt(2) part.
B[i,j]:=w * Eltseq(T[i,j])[2] * w; //Eltseq*w is the number. Then multiplying it by w again makes it rational.
end if; //b/c apparently just multiplying it by w is enough to give the rational part.
end for;
end for;
B:=ChangeRing(B, Integers());
return B;
end function;

Nebegram:=function(m); //The gram matrix
G:=Nebegen(m);
G:=G*Transpose(G);
return G;
end function;

Griessgen:=function(m); //Computes the generator matrix for Griess's constrction by rescaling Nebe et al's matrix.
Q<w>:=QuadraticField(2);
g:=Nebegen(m);
g:=ChangeRing(g,QuadraticField(2));
if IsEven(m) then
c:=w ^ (-m div 2);
end if;
if IsOdd(m) then
c:=w ^ ((-m-1) div 2);
end if;
for i:=1 to 2^m do
for j:=1 to 2^m do
g[i,j]:=c * g[i,j];
end for;
end for;
return g;
end function;

Griessgram:=function(m); //This one calculates the Gram matrix using the results of the previous function
curlym:=Griessgen(m);
curlym:=curlym*Transpose(curlym);
curlym:=ChangeRing(curlym,Integers());
return curlym;
end function;

UnimodularBRW:=function(m); //This gives a unimodular representation of BRW by conjugating each generator of BRW by the generator matrix of the Barnes-Wall lattice.
N:=ChangeRing(Nebegen(m),RationalField());
L:=[];
for A in Generators(BRW(m)) do
B:=ChangeRing(A,RationalField());
L:=Append(L,ChangeRing(N*B*N^(-1),Integers()));
end for;
return sub<GeneralLinearGroup(2^m,Integers())|L>;
end function;

PrimeBRW:=function(m,p);
M:=[];
for x in Generators(UnimodularBRW(m)) do
M:=Append(M,ChangeRing(x,FiniteField(p)));
end for;
return sub<GeneralLinearGroup(2^m,FiniteField(p)) | M>;
end function;

Grassmannian:=function(d,n); //This computes the number of d-dimensional subspaces of an n-dimensional vector space over Z2.
a:=2^n-1;
b:=2^d-1;
for i:=1 to d-1 do
a:=a*(2^n-2^i);
b:=b*(2^d-2^i);
end for;
return a/b;
end function;

Centralizer2:=function(m); //This computes the centralizer of a non-identity element of order 2 element in E(m). They're all isomorphic so I just get one.
G:=E(m);
Two:=0;
for x in G do
if Order(x) eq 2 and not x in Center(G) then
Two:=x;
break;
end if;
end for;
return Centralizer(G,Two);
end function;

Centralizer4:=function(m); //Same as above, for order 4.
G:=E(m);
Four:=0;
for x in G do
if Order(x) eq 4 then
Four:=x;
break;
end if;
end for;
return Centralizer(G,Four);
end function;