balanced:=function(m);
Q<w>:=QuadraticField(2);
M:=Matrix([[w,0],[1,1]]);
T:=M;
for i:=1 to m-1 do// This performs the tensor product m times.
T:=TensorProduct(T,M);
end for;
return T;
end function;

Nebegen:=function(m);
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);
G:=Nebegen(m);
G:=G*Transpose(G);
return G;
end function;

Griessgen:=function(m);
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;

Griessgramdirect:=function(m); //This one calculates it directly.
G:=ScalarMatrix(2 ^ m,0); //Initializing the Gram matrix.
for a:=1 to 2^m do //Looping through the entries of G.
for b:=1 to 2^m do
i:=a - 1; //I have to do this b/c our formula has indices start at 0.
j:=b - 1;
s:=0; //This will become the i,jth entry of the Gram matrix.
x:=0; //This will be the number of times the binary reps of i and j match in a slot.
z:=0; //This will be the number of zeroes in the binary rep of i.
Z:=0; //This will be the number of zeroes in the binary rep of j.
t:=0; //This will me which equation to use: 0 if i,j even; 1 if i or j odd; 2 if i,j odd.
L:=[* *]; //Initializing the binary rep of i.
M:=[* *]; //Initializing the binary rep of j.
for l:=1 to m do
if (i - (2 ^ (m - l))) ge 0 then
i:=i - (2 ^ (m - l));
L:=Append(L,1);
else
L:=Append(L,0);
end if;
if (j - (2 ^ (m - l))) ge 0 then
j:=j - (2 ^ (m - l));
M:=Append(M,1);
else
M:=Append(M,0);
end if;
end for;
for l:=1 to (m) do
if L[l] eq 0 then
z:=z+1;
end if;
if M[l] eq 0 then
Z:=Z+1;
end if;
if L[l] eq M[l] then
x:=x+1;
end if;
end for;
if IsOdd(z) then
t:=t+1;
end if;
if IsOdd(Z) then
t:=t+1;
end if;
if t eq 0 then
s:= 2 ^ (x div 2);
elif t eq 1 then
s:= 2 ^ ((x + 1) div 2);
elif t eq 2 then
s:= 2 ^ ((x div 2) +1 );
end if;
G[a,b]:=s;
end for;
end for;
return G;
end function;

//This function sends a complex number to a 2x2 real matrix.
Realify:=function(j);
Q<i>:=QuadraticField(-1);
j:=j + (0 * i);
M:=Matrix([[Eltseq(j)[1],Eltseq(j)[2]],[-1 * Eltseq(j)[2],Eltseq(j)[1]]]);
M:=ChangeRing(M, Integers());
return M;
end function;

//This function sends a complex n x n matrix to a 2n x 2n real matrix.
RealifyMatrix:=function(M);
Q<i>:=QuadraticField(-1);
M:=ChangeRing(M,QuadraticField(-1));
R:=ScalarMatrix(2 * NumberOfRows(M), 0); //Initializing the real matrix.
l:=-1;
k:=-1;
for a:=1 to NumberOfRows(M) do //Looping through the rows of the complex matrix
L:=l + (2 * a);
for b:=1 to NumberOfRows(M) do
K:=k + (2 * b);
R:=InsertBlock(R,Realify(M[a,b]),L,K);
end for;
end for;
return R;
end function;

//This function gets the IEEE construction of the Barnes Wall lattice of degree m, via the definition; start with {{1,1},{0,1+i}} and tensor m times.
IEEEimaginary:=function(m);
Q<i>:=QuadraticField(-1);
M:=Matrix([[1,1],[0,1+i]]);
T:=M;
for i:=1 to m-1 do
T:=TensorProduct(T,M);
end for;
return T;
end function;

//This function realifies the complex generator matrix given by the IEEE construction.
IEEEreal:=function(m);
Q<i>:=QuadraticField(-1);
R:=ChangeRing(RealifyMatrix(IEEEimaginary(m)),Integers());
return R;
end function;

//This function gets the MathOverflow construction of the Barnes Wall lattice of degree m, via the definition; starting with {{1,1},{1,i}} and tensoring m times.
Overflowimaginary:=function(m);
Q<i>:=QuadraticField(-1);
M:=Matrix([[1,1],[1,i]]);
T:=M;
for i:=1 to m-1 do
T:=TensorProduct(T,M);
end for;
return T;
end function;

//This function realifies the complex generator matrix given by the Overflow construction.
Overflowreal:=function(m);
Q<i>:=QuadraticField(-1);
R:=ChangeRing(RealifyMatrix(Overflowimaginary(m)),Integers());
return R;
end function;


//This gets the Gram matrix of the complex generator matrix given by the IEEE construction.
IEEEimaginarygram:=function(m);
X:=IEEEimaginary(m);
X:=X*Transpose(X);
return X;
end function;

//This gets the Gram matrix of the realified generator matrix given by the IEEE construction.
IEEErealgram:=function(m);
X:=IEEEreal(m);
X:=X*Transpose(X);
return X;
end function;

//The gram matrix of the IEEE realified matrix is NOT the realified version of the complex gram matrix; this will return false
RealifyMatrix(IEEEimaginarygram(1)) eq IEEErealgram(1);

//This gets the Gram matrix of the realified generator matrix given by the Overflow construction.
Overflowrealgram:=function(m);
X:=Overflowreal(m);
X:=X*Transpose(X);
return X;
end function;

//This gets the Gram matrix of the complex generator matrix given by the Overflow construction.
Overflowimaginarygram:=function(m);
X:=Overflowimaginary(m);
X:=X*Transpose(X);
return X;
end function;

//This directly calculates the complex generator matrix of IEEE's construction. If ever ik=1 and jk=0, that entry is 0. Otherwise, the ij'th entry is (1+i) to the power of the number of times ik=jk=1.
Ione:=function(m);
Q<i>:=QuadraticField(-1);
G:=ScalarMatrix(2^m,i);
for a:=1 to 2^m do
for b:=1 to 2^m do
z:=a - 1; //z instead of i since i is in use
j:=b - 1;
s:=1;
x:=0;
L:=[* *];
M:=[* *];
for l:=1 to m do
if (z - (2 ^ (m - l))) ge 0 then
z:=z - (2 ^ (m - l));
L:=Append(L,1);
else
L:=Append(L,0);
end if;
if (j - (2 ^ (m - l))) ge 0 then
j:=j - (2 ^ (m - l));
M:=Append(M,1);
else
M:=Append(M,0);
end if;
end for;
for l:=1 to m do
if L[l] eq 1 then
if M[l] eq 0 then
s:=0;
elif M[l] eq 1 then
x:=x+1;
end if;
end if;
end for;
if s ne 0 then
s:=(1 + i) ^ x;
end if;
G[a,b]:=s;
end for;
end for;
return G;
end function;

//This directly calculates the complex generator matrix of MathOverflow's construction. The ij'th entry is i to the power of the number of times ik=jk=1.
Itwo:=function(m);
Q<i>:=QuadraticField(-1);
G:=ScalarMatrix(2^(m),i);
for a:=1 to 2^m do
for b:=1 to 2^m do
z:=a - 1; //z instead of i since i is in use
j:=b - 1;
x:=0;
L:=[* *];
M:=[* *];
for l:=1 to m do
if (z - (2 ^ (m - l))) ge 0 then
z:=z - (2 ^ (m - l));
L:=Append(L,1);
else
L:=Append(L,0);
end if;
if (j - (2 ^ (m - l))) ge 0 then
j:=j - (2 ^ (m - l));
M:=Append(M,1);
else
M:=Append(M,0);
end if;
end for;
for l:=1 to m do
if L[l] eq 1 then
if M[l] eq 1 then
x:=x+1;
end if;
end if;
end for;
G[a,b]:=i ^ x;
end for;
end for;
return G;
end function;

//This directly calculates the realifed generator matrix of MathOverflow's construction. Take i^x to the power of the number of times that the binary reps of i and j coincide in all but the last slot, realify it, and then pick the entry of that matrix corresponding to the last slots of the binary reps of i and j.
Rtwo:=function(k);
m:=k+1;
Q<i>:=QuadraticField(-1);
G:=ScalarMatrix(2^(m),i);
for a:=1 to 2^(m) do
for b:=1 to 2^(m) do
z:=a - 1; //z instead of i since i is in use
j:=b - 1;
s:=0;
x:=0;
L:=[* *];
M:=[* *];
for l:=1 to m do
if (z - (2 ^ (m - l))) ge 0 then
z:=z - (2 ^ (m - l));
L:=Append(L,1);
else
L:=Append(L,0);
end if;
if (j - (2 ^ (m - l))) ge 0 then
j:=j - (2 ^ (m - l));
M:=Append(M,1);
else
M:=Append(M,0);
end if;
end for;
for l:=1 to k do
if L[l] eq 1 then
if M[l] eq 1 then
x:=x+1;
end if;
end if;
end for;
s:=i ^ x;
S:=Realify(s);
G[a,b]:=S[L[m]+1,M[m]+1];
end for;
end for;
G:=ChangeRing(G,Integers());
return G;
end function;

//This directly calculates the realifed generator matrix of IEEE's construction, using a similar method as before.
Rone:=function(k);
m:=k+1;
Q<i>:=QuadraticField(-1);
G:=ScalarMatrix(2^m,i);
for a:=1 to 2^m do
for b:=1 to 2^m do
z:=a - 1; //z instead of i since i is in use
j:=b - 1;
s:=1;
x:=0;
L:=[* *];
M:=[* *];
for l:=1 to m do
if (z - (2 ^ (m - l))) ge 0 then
z:=z - (2 ^ (m - l));
L:=Append(L,1);
else
L:=Append(L,0);
end if;
if (j - (2 ^ (m - l))) ge 0 then
j:=j - (2 ^ (m - l));
M:=Append(M,1);
else
M:=Append(M,0);
end if;
end for;
for l:=1 to k do
if L[l] eq 1 then
if M[l] eq 0 then
s:=0;
elif M[l] eq 1 then
x:=x+1;
end if;
end if;
end for;
if s ne 0 then
s:=(1 + i) ^ x;
end if;
S:=Realify(s);
G[a,b]:=S[L[m]+1,M[m]+1];
end for;
end for;
G:=ChangeRing(G,Integers());
return G;
end function;

//This function checks if the formulas used to directly calculate the matrices (looking at binary reps) match the matrices given by the definition (tensoring then realifying). All of them will be true for any m.
checkformulaswork:=function(m);
Ione(m) eq IEEEimaginary(m);
Itwo(m) eq Overflowimaginary(m);
Rone(m) eq IEEEreal(m);
Rtwo(m) eq Overflowreal(m);
return "";
end function;

//This function checks if the lattices generated by the realified generator matrices of the IEEE construction and Overflow construction match. It will return true for any m.
checklatticessame:=function(m);
return Lattice(Rone(m)) eq Lattice(Rtwo(m));
end function;

//This directly calculates the matrix that relates the complex generating matrices given by IEEE and Overflow's constructions.
ImaginaryRelatingMatrix:=function(m);
return Ione(m)*Itwo(m)^-1;
end function;

//This calculates the matrix that relates them via starting with the matrix for m=1, and tensoring m times.
ImaginaryRelatingMatrixTensor:=function(m);
M:=ImaginaryRelatingMatrix(1);
T:=M;
for i:=1 to m-1 do
T:=TensorProduct(T,M);
end for;
return T;
end function;

//This verifies that the above two functions give the same matrix
ImaginaryRelatingMatricesSame:=function(m);
return ImaginaryRelatingMatrix(m) eq ImaginaryRelatingMatrixTensor(m);
end function;

//This realifies the matrix that relates the complex generating matrix. As of this moment in time, we have no idea if the output of this is the matrix that relates the realified generating matrices, but it will turn out to be.
RealifyImaginaryRelatingMatrix:=function(m);
return RealifyMatrix(ImaginaryRelatingMatrixTensor(m));
end function;

//This directly calculates the matrix that relates the realified generating matrices.
RealRelatingMatrix:=function(m);
return Rone(m)*ChangeRing(Rtwo(m),RationalField())^-1;
end function;

//This checks to see if the matrices given by the above two functions are the same.
RelatingMatricesSame:=function(m)
return RealRelatingMatrix(m) eq RealifyImaginaryRelatingMatrix(m);
end function;

//For the IEEE and Overflow realified generating matrices to generate the same Z-lattice, their relating matrix has to be unimodular. This checks that. It will always return Integer Ring,1.
IsRelatingMatrixUnimodular:=function(m);
I:=RealifyImaginaryRelatingMatrix(1);
return BaseRing(I),Determinant(I);
end function;

//Basically what I've determined here is that IEEE and Overflow have the same realified Gram matrices, and they do generate the same realified lattice; the unimodular matrix that transforms one into the other is {{1,0},{i,-i}} tensored with itself m times and then realified.

SplitMatrix:=function(M);
Q<i>:=QuadraticField(-1);
R:=ScalarMatrix(NumberOfRows(M), 0);
I:=ScalarMatrix(NumberOfRows(M), 0);
for i:=1 to NumberOfRows(M) do
for j:=1 to NumberOfColumns(M) do
R[i,j]:=Eltseq(M[i,j])[1];
I[i,j]:=Eltseq(M[i,j])[2];
end for;
end for;
return <R,I>;
end function;

OtherRealifyMatrix:=function(M);
Q<i>:=QuadraticField(-1);
L:=ScalarMatrix(2 * NumberOfRows(M),0);
R:=SplitMatrix(M)[1];
I:=SplitMatrix(M)[2];
L:=InsertBlock(L,R,1,1);
L:=InsertBlock(L,I,1,NumberOfRows(M)+1);
L:=InsertBlock(L,-1*I,NumberOfColumns(M)+1,1);
L:=InsertBlock(L,R,NumberOfRows(M)+1,NumberOfRows(M)+1);
return L;
end function;

IEEEconstant:=function(m);
Q<w>:=QuadraticField(2);
if IsOdd(m) then
l:=(m-1) / 4;
else
l:=(m-2) / 4;
end if;
if IsCoercible(Integers(),l) then
c:=2^(-(Integers() ! l));
else
k:=Integers() ! (l-1/2);
c:=2^(-k) * 1/w;
end if;
return c;
end function;

CheckIsometric:=function(m);
G:=Griessgram(m);
if m eq 1 then
I:=Matrix([[1,0],[0,1]]);
else
I:=IEEEreal(m-1);
end if;
c:=IEEEconstant(m);
J:=c*I;
H:=J*Transpose(J);
H:=ChangeRing(H,Integers());
return IsIsometric(LatticeWithGram(G),LatticeWithGram(H));
end function;