/* V is the 4-dim vector space over GF(3) with alternating pairing
given by
<e1,e2> = 1, and <e3, e4> = 1 */

F := Matrix(GF(3),4,4,[[0,1,0,0],[-1,0,0,0],[0,0,0,1],[0,0,-1,0]]);
V := VectorSpace(GF(3),4,F);

/* G is the symplectic group */

G := IsometryGroup(V);


Lines := {@ {x,-x} : x in V | x ne Zero(V) @};

/*  Finding the set of all pairs {W, W^perp}, where
W is a non-isotropic plane  */

Planes := {@ sub<V|x,y> : x, y in V | DotProduct(x,y) ne 0 @};
Orderedplanepairs := {@ @};
Planepairs := {@ @};
for i := 1 to #Planes do
    W := Planes[i];
    Wperp := OrthogonalComplement(V,W);
    bool := true;
    for j := 1 to i-1 do
	if Planes[j] eq Wperp then
	    bool := false;
	    break;
	end if;
    end for;

    if bool then
	Include(~Orderedplanepairs,{@ W, Wperp @});
	Include(~Planepairs,{W,Wperp});
    end if;
end for;


/* Constructing the 85 x 130 intersection matrix defined in
Theorem 4.9 of "The Siegel Modular Variety of Degree 2 and level 3"
by Hoffman and Weintraub  */

temp_list := [];
for i := 1 to 85 do
    temp := [];
    if i le 45 then
	for j := 1 to 130 do
	    if j le 90 then
		if j mod 45 eq i mod 45 then
		    Append(~temp,-1);
		else
		    Append(~temp,0);
		end if;
	    else
		Append(~temp,0);
	    end if;
	end for;
    else
	line := Lines[i-45];
	for j := 1 to 130 do
	    if j le 45 then
		Delta := Orderedplanepairs[j];
		Wperp := Delta[2];
		if line subset Wperp then
		    Append(~temp,1);
		else
		    Append(~temp,0);
		end if;
	    elif j le 90 then
		Delta := Orderedplanepairs[j-45];
		W := Delta[1];
		if line subset W then
		    Append(~temp,1);
		else
		    Append(~temp,0);
		end if;
	    else
		second_line := Lines[j-90];
		if line eq second_line then
		    Append(~temp,-2);
		elif DotProduct(Random(line),Random(second_line)) eq 0 then
		    Append(~temp,1);
		else
		    Append(~temp,0);
		end if;
	    end if;
	end for;
    end if;
    Append(~temp_list,temp);
end for;
Intersection_matrix := Matrix(IntegerRing(),85,130,temp_list);


/* Defining the permutation matrices for the action of the two
generators of G on the free Z-module of dim 85 generated by Lines
and Planepairs {W,W^perp} */

Perm_mats := [];
for i := 1 to #Generators(G) do
    g := G.i;

    temp_list := [];

    for j := 1 to #Planepairs do
	Planepair := Planepairs[j];
	x := Random(Planepair);
	xi := x*Transpose(g);     /* We need to look
at the left regular reps Z[G/G_{40}] and Z[G/G_{45}] */
	n := Index(Planepairs,{xi,OrthogonalComplement(V,xi)});
	temp := [];
	for k := 1 to 85 do
	    if k ne n then
		Append(~temp,0);
	    else
		Append(~temp,1);
	    end if;
	end for;
	Append(~temp_list,temp);
    end for;

    for j := 1 to #Lines do
	Line := Lines[j];
	x := Random(Line);
	xi := x*Transpose(g);
	n := Index(Lines,{xi,-xi});
	temp := [];
	for k := 1 to 85 do
	    if k ne 45+n then
		Append(~temp,0);
	    else
		Append(~temp,1);
	    end if;
	end for;
	Append(~temp_list,temp);
    end for;

    perm_mat := Matrix(IntegerRing(),85,85,temp_list);
    Append(~Perm_mats,perm_mat);
end for;

/* FreeM is the free Z-module of dim 85 generated by Lines
and Planepairs */

FreeM := Domain(Intersection_matrix);
kernelsub := Kernel(Intersection_matrix);

/* Checking that the kernel is G-invariant  */

bool := true;
for g in Perm_mats do
    for i := 1 to Dimension(kernelsub) do
	if not kernelsub.i*g in kernelsub then
	    bool := false;
	    print "Bad";
	end if;
    end for;
end for;


/* Obtaining the quotient module M of dimension 61, and finding
matrices describing action of G on M. The module M is
H^2(A_2(3)*,Z) as per Proposition 1.4 of Hoffman and Weintraub */

M, q := quo<FreeM | kernelsub>;

Desired_rep_mats := [];
for i := 1 to #Generators(G) do
    g := G.i;

    temp_list := [];

    for j := 1 to Dimension(M) do
	x := M.j;
	xi := q((x @@ q)*Perm_mats[i]);
	Append(~temp_list,xi);
    end for;

    perm_mat := Matrix(IntegerRing(),61,61,temp_list);
    Append(~Desired_rep_mats,perm_mat);
end for;

x := Desired_rep_mats[1];
y := Desired_rep_mats[2];

/* OK, now one defines the group G in GL_61(Z) in terms of the
representation PSp_4(F_3) = <x,y> as above. */

G<x,y>:=MatrixGroup<61,Integers()|x,y>;

/* the explit matrices are as follows: */

print("##########################");
print("begin first computation");
print("the first part of this code computes the action of G=PSp_4(F_3)");
print("on M = H^2(X,Z) \simeq Z^61. The explicit matrices x and y giving");
print("this action are printed below");
print("[Warning! matrices in Magma act on the right!]");
x;
y;


/* And finally, here is H^2(A(3)^*,Z) as a G-module! */

M:=GModule(G); MD:=Dual(M);

/* We already know what it is rationally, so just checking that is
correct */

X:=CharacterTable(G);

print("Here are some representations of G");
print("[Q]=[V_1]");
X[1];
print("[V_15]");
X[8];
print("[V_20]");
X[9];
print("[chi_24]");
X[10];
print("The permutation representation Q[G/G40]=[Q]+[chi_24]+[V_15]");
X[1]+X[10]+X[8];
print("The permutation representation Q[G/G45]=[Q]+[chi_24]+[V_20]");
X[1]+X[10]+X[9];
print("The representation Q[G/G45] + Q[G/G40] - [chi_24]");
(X[1]+X[10]+X[8])+(X[1]+X[10]+X[9]) - X[10];
print("This should agree with [M tensor Q] which is");
Character(M);

/* first computation, determine for which subgroups P of G one has
H^1(P,M) =/= 1 */

print("##########################");
print("In this computation, we determine for which subgroups P of G");
print("One has H^1(P,M) =/= 1");
print("##########################");

CM:=CohomologyModule(G,M);
CMD:=CohomologyModule(G,MD);

SetSeed(1);
SL:=SubgroupLattice(G);
ZG:=Subgroups(G);


for n in [2..116] do
P:=ZG[n]`subgroup;
CP:=Restriction(CM,P);
CPD:=Restriction(CMD,P);
count:=0;
if #CohomologyGroup(CP,1) ge 2 then
	count:=count+1;
	end if;
if #CohomologyGroup(CPD,1) ge 2 then
        count:=count+1;
        end if;
if count ge 1 then
	print("********************");
	print("for the subgroup H of label n, with [n,|H|,[G;H]] equal to");
	[n,Order(P),Index(G,P)];
	print("the group H^1(H,M) is equal to");
	CohomologyGroup(CP,1);
        print("the group H^1(H,M^vee) is equal to");
        CohomologyGroup(CPD,1);
end if;
end for;


print("##########################");
print("In this second computation, for all H in G, find lcm of H^1(P,M) for");
print("P inside H, and lcm of H^1(P,M^vee) for P inside H");
print("##########################");

/* second computation. For all H in G, find lcm of #H^1(P,M) for P inside H  */

for n in [2..115] do
	H:=ZG[n]`subgroup;
	print("*******************");
	print("for the subgroup H of label n, with [n,|H|,[G;H]] equal to");
        [n,Order(H),Index(G,H)];
	IdentifyGroup(H);
	L:=1; LD:=1;
	ZH:=Subgroups(H);
	for m in [2..#ZH] do
		P:=ZH[m]`subgroup;
		CP:=Restriction(CM,P);
		CPD:=Restriction(CMD,P);
		L:=LCM(L,GCD(6,#CohomologyGroup(CP,1)));
		LD:=LCM(LD,GCD(6,#CohomologyGroup(CPD,1)));
		end for;
		[L,LD];
	end for;


/* list the lattice of subgroups */

print("##########################");
print("list the lattice of subgroups of G = PSp_4(F_3)");
print("##########################");


SL;

print("##########################");
print("For every H in G, compute the Burnside Cokernel which is");
print("R_H(G) modulo permutation representations. Then determine");
print("if V_H where V=[M tensor Q] is non-trivial in this group,");
print("equivalently, determine if V_H is a virtual sum of permutation");
print("representations of H. If not, then print a message");
print("##########################");


X:=CharacterTable(G);
V:=X[1] + X[9] + X[10] + X[1] + X[8];
print("V has dimension");
V[1];
print("G has Burnside Cokernel");
BurnsideCokernel(G);


for m in [2..#ZG] do
n:=#ZG+2-m;
counter:=0;
H:=ZG[n]`subgroup;
VH:=Restriction(V,H);
if Index(G,H) ne 1080 then
	if IsPermutationCharacter(VH) then
	counter:=counter+1;
	else
	print("************************");
	print("For the subgroup H of label n then [n,[G:H]] equals");
	[n,Index(G,H)];
	if Order(H) le 1000 then
		IdentifyGroup(H);
	end if;
	print("X_H is not rational if the following is false");
	IsPermutationCharacter(VH);
	print("the order of V_H in the Burnside cokernel is as follows");
	counter:=0;
	for k in [1..6] do
		if counter eq 0 then
			if IsPermutationCharacter(k*VH) then
				counter:=k;
				end if;
			end if;
		end for;
	counter;
end if;
end if;
end for;