/*

The command FindRankFq(p,n,m,r,k) returns a sequece consists of the Q-rank, Fq-rank where q<=p and q is a prime for the group B4(Zn+Zm+Zr+Zk).


*/






// Create the function checking if (i1,j1,k1,l1),(i2,j2,k2,l2),(i3,j3,k3,l3),(i4,j4,k4,l4) generates Zn+Zm+Zr+Zk.
CheckGenerators:=function (n,m,r,k,i1,j1,k1,l1,i2,j2,k2,l2,ix3,j3,k3,l3,i4,j4,k4,l4);
if [n,m,r] eq [1,1,1] then//Cyclic groups
if Gcd(Gcd(Gcd(Gcd(k,l1),l2),l3), l4) eq 1 then
return true;
else
return false;
end if;
elif [n,m] eq [1,1] then

Z:=RModule(IntegerRing(),2);
L:=sub<Z|[r,0],[0,k]>;
Znm:=quo<Z|L>;
gen:=sub<Znm|[k1,l1],[k2,l2],[k3,l3],[k4,l4]>;

if gen eq Znm then
return true;
else
return false;
end if;

elif n eq 1 then

Z:=RModule(IntegerRing(),3);
L:=sub<Z|[m,0,0],[0,r,0],[0,0,k]>;
Znm:=quo<Z|L>;
gen:=sub<Znm|[j1,k1,l1],[j2,k2,l2],[j3,k3,l3],[j4,k4,l4]>;

if gen eq Znm then
return true;
else
return false;
end if;

else
Z:=RModule(IntegerRing(),4);
L:=sub<Z|[n,0,0,0],[0,m,0,0],[0,0,r,0],[0,0,0,k]>;
M:=quo<Z|L>;
gen:=sub<M|[i1,j1,k1,l1],[i2,j2,k2,l2],[i3,j3,k3,l3],[i4,j4,k4,l4]>;

if gen eq M then
return true;
else
return false;
end if;
end if;
end function;




// Turn Zn+Zm+Zr+Zk to a standard form
StandardForm:=function(n,m,r,k);
S:=Sort([n,m,r,k]);
n:=S[1];
m:=S[2];
r:=S[3];
k:=S[4];
prod:=n*m;
n:=Gcd(m,n);
m:=prod div n;
prod:=m*r;
m:=Gcd(m,r);
r:=prod div m;
prod:=r*k;
r:=Gcd(r,k);
k:=prod div r;

return Sort([n,m,r,k]);

end function;






//Compute the matrix and the rank
FindRankFq:=function(q,n,m,r,k);

//Turn n,m,r,k into standard forms
Order:=StandardForm(n,m,r,k);
n:=Order[1];
m:=Order[2];
r:=Order[3];
k:=Order[4];

//Store Generators of B4(Zn+Zm+Zr+Zl)
B4Gens:=AssociativeArray();
B4NGens:=0; //number of generators 
B4Index:=AssociativeArray(RModule(IntegerRing(),16)); //record the order of generators

for i1,i2,i3,i4 in [0..n-1] do
for j1,j2,j3,j4 in [0..m-1] do
for k1,k2,k3,k4 in [0..r-1] do
for l1,l2,l3,l4 in [0..k-1] do
if CheckGenerators(n,m,r,k,i1,j1,k1,l1,i2,j2,k2,l2,i3,j3,k3,l3,i4,j4,k4,l4) then
B4NGens:=B4NGens+1;
B4Gens[B4NGens]:=[i1,j1,k1,l1,i2,j2,k2,l2,i3,j3,k3,l3,i4,j4,k4,l4];
B4Index[[i1,j1,k1,l1,i2,j2,k2,l2,i3,j3,k3,l3,i4,j4,k4,l4]]:=B4NGens;
end if;
end for;
end for;
end for;
end for;

row_number:=B4NGens*25;
currentrow:=0;
a:=SparseMatrix(RationalField(), row_number, B4NGens);
for i in [1..B4NGens] do
//Symmetry relation (S)
if (B4Gens[i][1..4] le B4Gens[i][5..8]) and (B4Gens[i][5..8] le B4Gens[i][9..12]) and (B4Gens[i][9..12] le B4Gens[i][13..16]) then//Find (a1,a2,a3,a4) satisfying lexiographic order

if (B4Gens[i][1..4] eq B4Gens[i][13..16]) then//No need for symmetry on (a,a,a,a)

elif (B4Gens[i][1..4] eq B4Gens[i][9..12]) and (B4Gens[i][9..12] lt B4Gens[i][13..16])then//Cases when (a,a,a,b)=(a,a,b,a)=(a,b,a,a)=(b,a,a,a) a<b
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][1..4] cat B4Gens[i][1..4] cat B4Gens[i][13..16] cat B4Gens[i][1..4];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;

currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][1..4] cat B4Gens[i][13..16] cat B4Gens[i][1..4]  cat B4Gens[i][1..4];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;

currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][13..16] cat B4Gens[i][1..4] cat B4Gens[i][1..4] cat B4Gens[i][1..4];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;


elif (B4Gens[i][1..4] lt B4Gens[i][5..8]) and (B4Gens[i][5..8] eq B4Gens[i][13..16])then//Cases when (b,a,a,a)=(a,a,b,a)=(a,b,a,a)=(a,a,a,b) a>b
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][5..8] cat B4Gens[i][5..8] cat B4Gens[i][1..4] cat B4Gens[i][5..8];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;

currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][5..8] cat B4Gens[i][1..4] cat B4Gens[i][5..8]  cat B4Gens[i][5..8];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;

currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][1..4] cat B4Gens[i][5..8] cat B4Gens[i][5..8] cat B4Gens[i][5..8];
index:=B4Index[target];
a[currentrow, index]:=a[currentrow, index]-1;


else // Other Cases 
for x in Sym(4) do 
if x ne Id(Sym(4)) then
indice1:=4*([1,2,3,4]^x)[1]-3;
indice2:=4*([1,2,3,4]^x)[2]-3;
indice3:=4*([1,2,3,4]^x)[3]-3;
indice4:=4*([1,2,3,4]^x)[4]-3;
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][indice1..indice1+3] cat B4Gens[i][indice2..indice2+3] cat B4Gens[i][indice3..indice3+3] cat B4Gens[i][indice4..indice4+3];
index:=B4Index[target];
a[currentrow,index]:=a[currentrow,index]-1;
end if;
end for;
end if;
end if;


//Relation (B)
if B4Gens[i][1..4] eq B4Gens[i][5..8] then
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
target:=B4Gens[i][1..4] cat [0,0,0,0] cat B4Gens[i][9..12] cat B4Gens[i][13..16];
index:=B4Index[target];
a[currentrow,index]:=a[currentrow,index]-1;
else
currentrow:=currentrow+1;
a[currentrow,i]:=a[currentrow,i]+1;
A1:=(B4Gens[i][1]-B4Gens[i][5]) mod n;
A2:=(B4Gens[i][2]-B4Gens[i][6]) mod m;
A3:=(B4Gens[i][3]-B4Gens[i][7]) mod r;
A4:=(B4Gens[i][4]-B4Gens[i][8]) mod k;
A:=[A1, A2, A3, A4];
target1:=A cat B4Gens[i][5..8] cat B4Gens[i][9..12] cat B4Gens[i][13..16];
B:=[(-A1) mod n, (-A2) mod m, (-A3) mod r, (-A4) mod k];
target2:=B4Gens[i][1..4] cat B cat B4Gens[i][9..12] cat B4Gens[i][13..16];
index1:=B4Index[target1];
index2:=B4Index[target2];
a[currentrow,index1]:=a[currentrow,index1]-1;
a[currentrow,index2]:=a[currentrow,index2]-1;
end if;
end for;


Result:=AssociativeArray();
Result[1]:=B4NGens-Rank(a);
primecount:=1;

for p in [1..q] do
if IsPrime(p) then
primecount:=primecount+1;
b:=ChangeRing(a,FiniteField(p)); 
Result[primecount]:=B4NGens-Rank(b);
end if;
end for;

result:=Matrix(RationalField(), 1, primecount, [<1,i, Result[i]>: i in [1..primecount]]);

return result;
end function;