Paramodularity for C277

GrpExt.m

+//////////////////////////////////////////////////
+// Implementation of parabolic group extensions
+// Author: Gonzalo Tornaría <tornaria@cmat.edu.uy>
+//////////////////////////////////////////////////
+
+declare type GrpExt;
+
+declare attributes GrpExt :
+  group, abelian, self, ambient, pi, gmodule, gmap, _gmodules;
+
+intrinsic ParabolicExtension(G :: GrpMat, A :: ModMatFld) -> GrpExt
+  { Construct the parabolic extension of G by A }
+  one := Matrix(One(G)); zero := Matrix(Zero(A));
+  require IsCompatible(one, zero)
+          : "G and A must have equal degree and base ring";
+  require &and [ g*a*g^-1 in A
+               : g in Generators(G), a in Generators(A)]
+          : "G must act on A by conjugation";
+  n := Degree(G);
+  R := BaseRing(G);
+  GL2n := GL(2*n, R);
+  // c : A --> ext
+  c := func< a | BlockMatrix([[one, a], [zero, one]]) >;
+  Cs := [c(r*a) : a in OrderedGenerators(A), r in Basis(R)];
+  // r : G --> ext
+  r := func< g | BlockMatrix([[g, zero], [zero, g]]) >;
+  Rs := [r(g) : g in OrderedGenerators(G)];
+  // p : ext --> G
+  p := func< x | ExtractBlock(x, 1, 1, n, n) >;
+  //
+  e := New(GrpExt);
+  e`group := G;
+  e`abelian := sub< GL2n | Cs >;
+  e`self := sub< GL2n | Cs cat Rs >;
+  e`ambient := e;
+  Ps := [p(x) : x in OrderedGenerators(e`self) ];
+  e`pi := hom< e`self -> e`group | Ps >;
+  cinv := func< x | ExtractBlock(x, 1, 1+n, n, n) >;
+  MA := MatrixAlgebra(R, Dimension(A));
+  gmap := func< x | Coordinates(A, A!cinv(x)) >;
+  gmat := func< g | MA ! [gmap(g^(-1)*c(a)*g) : a in Basis(A)] >;
+  action := [MA | gmat(r(g)) : g in Generators(G)];
+  // gmodule corresponding to the ambient
+  e`gmodule := GModule(G, action);
+  // e`gmap : e`abelian -> e`gmodule
+  e`gmap := hom< e`abelian -> e`gmodule | x :-> gmap(x) >;
+  // cache gmodules in the ambient for consistency
+  e`ambient`_gmodules := AssociativeArray();
+  e`ambient`_gmodules[G] := e`ambient`gmodule;
+
+  // it is required that G acts on A by conjugation
+  assert Kernel(e`pi) eq e`abelian;
+  assert Image(e`pi) eq e`group;
+  assert #e`abelian eq #A;
+  assert #e`self eq #e`group * #e`abelian;
+  return e;
+end intrinsic
+
+function sub(e, s)
+  assert s subset e`self;
+  f := New(GrpExt);
+  f`group := e`pi(s);
+  f`abelian := e`abelian meet s;
+  f`self := s;
+  f`ambient := e`ambient;
+  Ps := [e`pi(x) : x in OrderedGenerators(f`self)];
+  // compute gmodule
+  f`pi := hom< f`self -> f`group | Ps>;
+  b, gm := IsDefined(f`ambient`_gmodules, f`group);
+  if not b then
+    gm := Restriction(f`ambient`gmodule, f`group);
+    f`ambient`_gmodules[f`group] := gm;
+  end if;
+  f`gmodule := sub< gm |
+         [f`ambient`gmap(a) : a in Generators(f`abelian)] >;
+  f`gmap := hom< f`abelian -> f`gmodule | x :-> f`ambient`gmap(x) >;
+  assert Kernel(f`pi) eq f`abelian;
+  assert Image(f`pi) eq f`group;
+  assert #f`self eq #f`group * #f`abelian;
+  return f;
+end function;
+
+//////////////////////////////////////////////////
+
+// Some Lie algebras (gl_n, sl_n, sp_n, gsp_n)
+
+function gl(n, R)
+  return RMatrixSpace(R, n, n);
+end function;
+
+function sl(n, R)
+  A := gl(n, R);
+  B := RSpace(R, 1);
+  m := hom<A->B | [ [Trace(m)]
+                  : m in OrderedGenerators(A)]>;
+  return Kernel(m);
+end function;
+
+// the standard symplectic form used by Sp(n, q) in magma
+function symplectic_form(n, R)
+  assert IsEven(n);
+  J := Matrix(R, n, [i+j eq n+1
+                 select (i lt j select 1 else -1)
+                 else 0
+               : i,j in [1..n] ]);
+  assert J+Transpose(J) eq 0;
+  assert J*Transpose(J) eq 1;
+  return J;
+end function;
+
+function sp(n, R)
+  assert IsEven(n);
+  J := symplectic_form(n, R);
+  A := gl(n, R);
+  m := hom<A->A | [ mJ-Transpose(mJ)
+                    where mJ is m*J
+                  : m in OrderedGenerators(A)]>;
+  return Kernel(m);
+end function;
+
+function gsp(n, R)
+  assert IsEven(n);
+  J := symplectic_form(n, R);
+  assert J[1,n] eq 1;
+  A := gl(n, R);
+  m := hom<A->A | [ m1 - mu*J
+                    where mu := m1[1,n] // m[1,1] + m[n,n]
+                    where m1 := mJ-Transpose(mJ)
+                    where mJ is m*J
+                  : m in OrderedGenerators(A)]>;
+  return Kernel(m);
+end function;
+
+intrinsic SmallParabolic(G :: GrpMat) -> GrpExt
+  { Small parabolic as a group extension }
+  n := Degree(G);
+  R := BaseRing(G);
+  require IsEven(n) : "Degree of G must be even";
+  spn := sp(n, R);
+  gspn := gsp(n, R);
+  e := ParabolicExtension(G, spn);
+  // use the normalizer as ambient space
+  a := ParabolicExtension(G, gspn);
+  return sub(a, e`self);
+end intrinsic;
+
+intrinsic Normalizer(e :: GrpExt) -> GrpExt
+  { Normalizer of a group extension inside its ambient space }
+  return sub(e`ambient, Normalizer(e`ambient`self, e`self));
+end intrinsic;
+
+intrinsic Centralizer(e :: GrpExt) -> GrpExt
+  { Centralizer of Normalizer(e) }
+  return sub(e`ambient, Centralizer(Normalizer(e)`self, e`self));
+end intrinsic;
+
+intrinsic GModule(e :: GrpExt) -> ModGrp
+  { GModule corresponding to e }
+  return e`gmodule;
+end intrinsic;
+
+//////////////////////////////////////////////////
+
+intrinsic Group(e :: GrpExt) -> GrpMat
+  { .. }
+  return e`group;
+end intrinsic;
+
+intrinsic Abelian(e :: GrpExt) -> GrpMat
+  { .. }
+  return e`abelian;
+end intrinsic;
+
+intrinsic Self(e :: GrpExt) -> GrpMat
+  { .. }
+  return e`self;
+end intrinsic;
+
+intrinsic Ambient(e :: GrpExt) -> GrpMat
+  { .. }
+  return e`ambient;
+end intrinsic;
+
+intrinsic BaseRing(e :: GrpExt) -> Rng
+  { .. }
+  return BaseRing(Group(e));
+end intrinsic;
+
+intrinsic Dimension(e :: GrpExt) -> RngIntElt
+  { .. }
+  return Integers() ! Log(#BaseRing(e), #Abelian(e));
+end intrinsic;
+
+intrinsic '#'(e :: GrpExt) -> RngIntElt
+  { .. }
+  return #Self(e);
+end intrinsic;
+
+intrinsic GroupNameExt(G :: GrpMat : TeX := false) -> Str
+  { ... }
+  GN := GroupName(G : TeX:=TeX);
+  G_3 := [g : g in G | Order(g) eq 3];
+  G_6 := [g : g in G | Order(g) eq 6];
+  if (#G_3 gt 0 and &and [Trace(g) eq 0 : g in G_3])
+     or
+     (#G_6 gt 0 and &and [Trace(g) eq 0 : g in G_6])
+  then
+    return GN cat "(A)";
+  end if;
+  if (#G_3 gt 0 and &and [Trace(g) eq 1 : g in G_3])
+     or
+     (#G_6 gt 0 and &and [Trace(g) eq 1 : g in G_6])
+  then
+    return GN cat "(B)";
+  end if;
+  return GN;
+end intrinsic;
+
+intrinsic Print(e :: GrpExt)
+  { .. }
+  printf "GrpExt of %o by GF(%o)^%o",
+      GroupNameExt(Group(e)), #BaseRing(e), Dimension(e);
+end intrinsic;
+
+intrinsic Restriction(e :: GrpExt, H :: GrpMat) -> GrpExt
+  { .. }
+  require H subset e`group : "H must be a subgroup of G";
+  return sub(e, H @@ e`pi meet e`self);
+end intrinsic;
+
+intrinsic Conjugates(e :: GrpExt, f :: GrpExt) -> []
+  { .. }
+  return Conjugates(e`self, f`self);
+end intrinsic;
+
+intrinsic Conjugates(e :: GrpExt) -> []
+  { .. }
+  return Conjugates(e`ambient`self, e`self);
+end intrinsic;
+
+intrinsic IsConjugate(e :: GrpExt, f :: GrpExt) -> BoolElt
+  { .. }
+  require e`ambient eq f`ambient: "ambient spaces must be equal";
+  return IsConjugate(e`ambient`self, e`self, f`self);
+end intrinsic;
+
+intrinsic FilterConjugates(es :: [GrpExt]) -> []
+  { .. }
+  fs := [];
+  ii := [];
+  for i := 1 to #es do
+    for j := 1 to #fs do
+      if IsConjugate(es[i], fs[j]) then
+        Append(~ii[j], i);
+        continue i;
+      end if;
+    end for;
+    // not found
+    Append(~fs, es[i]);
+    Append(~ii, [i]);
+  end for;
+  return fs, ii;
+end intrinsic;
+
+intrinsic Subextensions(e :: GrpExt : Codim := -1) -> []
+  { .. }
+  if Codim eq -1 then
+    candidates := [ s`subgroup : s in
+                    Subgroups(e`self
+                             : OrderMultipleOf:=#e`group
+                             , OrderDividing:=#e`self) ];
+  elif Codim gt Dimension(e) then
+    return [];
+  else
+    q := #BaseRing(e);
+    o := IntegerRing() !  (#e/q^Codim);
+    candidates := [ s`subgroup : s in
+                    Subgroups(e`self : OrderEqual := o) ];
+  end if;
+  return [ sub(e, s) : s in candidates | e`pi(s) eq e`group ];
+end intrinsic;
+
+intrinsic 'meet' (e :: GrpExt, f :: GrpExt) -> GrpExt
+  { .. }
+  return sub(e, e`self meet f`self);
+end intrinsic;
+
+intrinsic 'subset' (e :: GrpExt, f :: GrpExt) -> BoolElt
+  { .. }
+  return Self(e) subset Self(f);
+end intrinsic;
+
+intrinsic 'eq' (e :: GrpExt, f :: GrpExt) -> BoolElt
+  { .. }
+  assert e`ambient`self eq f`ambient`self;
+  return Self(e) eq Self(f);
+end intrinsic;
+
+// m can be of type HomGrp or GrpAutoElt
+intrinsic IsIsomorphism(m, e :: GrpExt, f :: GrpExt) -> BoolElt
+  { .. }
+  return IsTrivial(Kernel(m)) and Image(m) eq f`self and
+        // m should commute with the projections pi
+        &and[ x @ e`pi eq x @ m @ f`pi : x in Generators(e`self) ];
+end intrinsic;
+
+intrinsic IsIsomorphic (e :: GrpExt, f :: GrpExt) -> BoolElt, Map
+  { test isomorphism as group extensions }
+  if e`group ne f`group or #e`abelian ne #f`abelian then
+    return false;
+  end if;
+
+  b, m := IsIsomorphic(e`self, f`self);
+
+  if not b then
+    return false;
+  end if;
+
+  // may or may not be isomorphic as group extensions
+  error if m(e`abelian) ne f`abelian, "Not implemented -- FIXME";
+
+  b, g := IsInnerAutomorphism(e`group, Inverse(e`pi) * m * f`pi);
+
+  // may or may not be isomorphic as group extensions
+  error if not b, "Not implemented for non-inner -- FIXME";
+
+  gg := g @@ f`pi;
+
+  mm := hom< e`self->f`self |
+             [gg*m(x)*gg^-1 : x in OrderedGenerators(e`self)] >;
+
+  // mm should be an isomorphism
+  assert IsIsomorphism(mm, e, f);
+
+  return true, mm;
+end intrinsic;
+
+intrinsic AutomorphismGroup(e :: GrpExt) -> GrpAuto
+  { automorphisms as group extensions }
+
+  aut := AutomorphismGroup(e`self);
+  autg := AutomorphismGroup(e`group);
+
+  // m1 : autfp -> aut
+  // m2 : autfp -> autp
+  autfp, m1 := FPGroup(aut);
+  autp, m2 := PermutationGroup(autfp);
+
+  k := Kernel(hom<autp->autg |
+     [ hom<e`group->e`group |
+          [x @@ e`pi @ (a@@m2@m1) @ e`pi : x in OrderedGenerators(e`group)]>
+       : a in OrderedGenerators(autp)] >);
+
+  return k;
+end intrinsic;
+
+intrinsic OuterAutomorphismRepresentatives(e :: GrpExt) -> []
+  { outer aut as gp ext (modulo inner aut by normalizer !!!)}
+
+  aut := AutomorphismGroup(e`self);
+  autg := AutomorphismGroup(e`group);
+
+  // m1 : autfp -> aut
+  // m2 : autfp -> outfp
+  // m3 : outfp -> outp
+  autfp, m1 := FPGroup(aut);
+  outfp, m2 := OuterFPGroup(aut);
+  outp, m3 := PermutationGroup(outfp);
+  // m0 : outp -> aut
+  m0 := Inverse(m3) * Inverse(m2) * m1;
+
+  out := [ m0(x) : x in outp ];
+
+  /* error if not &and[ IsIsomorphism(x, e, e) : x in out ], */
+  /*   "an outer automorphism is not identity on G -- FIXME"; */
+
+  return out;
+
+  /*
+  v := e`abelian;
+  vg := e`self;
+
+  out3 := [m0(x) : x in outp | v @ m0(x) eq v and vg @ m0(x) eq vg];
+
+  // must fix x as in 'IsIsomorphic' so it is an isomorphism
+  // however, it seems to be ok due to the way magma computes aut gp
+  out2 := [ m0(x) : x in outp | b where b,g := IsInnerAutomorphism(e`group,
+     hom<e`group->e`group |
+          [g @@ e`pi @ m0(x) @ e`pi : g in OrderedGenerators(e`group)]>
+          )];
+
+  print #Kernel(m2), #Set(out3), #out3, #outp, #out2;
+
+  error if not &and[ IsIsomorphism((x), e, e) : x in out2 ],
+    "an outer automorphism is not identity on G -- FIXME";
+
+  return out2;
+  //return "-", #out2, #out3, #outp, "-";
+  */
+
+end intrinsic;
+
+intrinsic GaloisExtension(R2 :: FldAb, VGs :: []) -> GrpExt
+  { HACK, compute
+    VGs := IsomorphismClasses(Subextensions(SmallParabolic(Image(rho)))); }
+
+  gal, r, s := AbsoluteGaloisGroup(R2);
+
+  for VG in VGs do
+    b, m := IsIsomorphic(gal, VG`self);
+    if b then
+      print Dimension(VG);
+    end if;
+  end for;
+
+  return 0;
+
+end intrinsic;
+
+intrinsic GaloisExtension(R2 :: FldAb, R1 :: FldNum, rho :: HomGrp) -> GrpExt
+  { Construct an extension corresponding to rho.
+      
+    WARNING: This is a hack that works by finding a parabolic extension so
+    it will be ok as long as G has no outer automorphisms, otherwise
+    we should fix it by checking traces of frobenius or something like
+    that.
+
+    IOW the aim at first is to construct an extension corresponding to
+    some conjugate of rho }
+
+  gal, r, s := AbsoluteGaloisGroup(R2);
+
+  G := Image(rho);
+  SG := SmallParabolic(G);
+  for VG in IsomorphismClasses(Subextensions(SG)) do
+    b, m := IsIsomorphic(gal, VG`self);
+    if b then
+      print Dimension(VG);
+    end if;
+  end for;
+
+  return 0;
+
+  /* one := Matrix(One(G)); zero := Matrix(Zero(A)); */
+  /* require Parent(one) eq Parent(zero) */
+  /*         : "G and A must have equal degree and base ring"; */
+  /* n := Degree(G); */
+  /* R := BaseRing(G); */
+  /* GL2n := GL(2*n, R); */
+  /* // c : A --> ext */
+  /* c := func< a | BlockMatrix([[one, a], [zero, one]]) >; */
+  /* Cs := [c(a) : a in OrderedGenerators(A)]; */
+  /* // r : G --> ext */
+  /* r := func< g | BlockMatrix([[g, zero], [zero, g]]) >; */
+  /* Rs := [r(g) : g in OrderedGenerators(G)]; */
+  /* // p : ext --> G */
+  /* p := func< x | ExtractBlock(x, 1, 1, n, n) >; */
+  /* // */
+  e := New(GrpExt);
+  e`group := Image(rho);
+  /* e`abelian := sub< GL2n | Cs >; */
+  /* e`self := sub< GL2n | Cs cat Rs >; */
+  /* e`ambient := e`self; */
+  /* Ps := [p(x) : x in OrderedGenerators(e`self) ]; */
+  /* e`pi := hom< e`self -> e`group | Ps >; */
+  /* // it is required that G acts on A by conjugation */
+  /* assert Kernel(e`pi) eq e`abelian; */
+  /* assert Image(e`pi) eq e`group; */
+  /* assert #e`self eq #e`group * #e`abelian; */
+  /* return e; */
+end intrinsic
+

Modularity.m

+////////////////////////
+// Proving modularity
+///////////////////////
+
+Attach("GrpExt.m");
+
+declare type Modularity;
+
+declare attributes Modularity :
+  R1, VGs, VHs, Vs, auts;
+
+// ASSUME
+//   i. rho : Gal(R1) -~-> G, a subgroup inside Sp(4,2)
+//  ii. G is abs irreducible in Sp(4,2)
+// iii. all automorphisms of G are inner [FIXME later]
+intrinsic ModNew (R1 :: FldNum, rho :: HomGrp) -> Modularity
+  { precomputation for modularity test }
+  G := Image(rho);
+  H := rho(Stabilizer(Domain(rho), 1));
+  SG := SmallParabolic(G);
+  VGs := FilterConjugates(Subextensions(SG));
+  auts := AssociativeArray();
+  VHs := AssociativeArray();
+  Vs := AssociativeArray();
+  // precompute automorphism groups
+  for VG in VGs do
+    auts[VG] := OuterAutomorphismRepresentatives(VG);
+    VHs[VG] := Restriction(VG, H);
+    Vs[VG] := VG`abelian;
+  end for;
+
+  m := New(Modularity);
+  m`R1 := R1;
+  m`VGs := VGs;
+  m`VHs := VHs;
+  m`Vs := Vs;
+  m`auts := auts;
+
+  return m;
+end intrinsic;
+
+intrinsic ModPhis(m :: Modularity, R2 :: FldAb) -> []
+  { return a list of all possible phi }
+
+  assert NumberField(BaseRing(R2)) eq m`R1;
+
+  // may not be the most efficient way to do this
+  //gal, r, s := AbsoluteGaloisGroup(R2);
+  gal, r, s := GaloisGroup(AbsoluteField(NumberField(R2)));
+
+  ans := [* *];
+
+  for VG in m`VGs do
+    b, phi := IsIsomorphic(gal, VG`self);
+    VH := m`VHs[VG]`self;
+    V := m`Vs[VG];
+    if b then
+      //print Dimension(VG);
+      for psi in m`auts[VG] do
+        assert psi(VH) eq VH;
+        assert psi(V) eq V;
+        Append(~ans, phi*psi);
+      end for;
+    end if;
+  end for;
+
+  return ans;
+
+end intrinsic;
+
+/////////////////////////
+// sorted degrees for a factorization (from high to low)
+degrees := func<fact | Reverse(Sort([Degree(f[1]): f in fact]))>;
+// return sorted degrees of factorization of pol mod p
+factpat := func<pol,p | degrees(Factorization(ChangeRing(pol,GF(p))))>;
+// turn a list of pairs <data,num> into a list repeating data num times
+untally := func<fs | [f[1] : j in [1..f[2]], f in fs]>;
+
+function tr(g)
+  return g[1,5]+g[2,6]+g[3,7]+g[4,8];
+end function;
+
+function tr1(g)
+  return g[1,1]+g[2,2]+g[3,3]+g[4,4];
+end function;
+
+ZZ := Integers();
+
+function show_classes(m,phi)
+  gal := Domain(phi);
+  cc := ConjugacyClasses(gal);
+  /**/
+  G := Group(m`VGs[1]);
+  cs := ConjugacyClasses(G);
+  cclass := func<mm | [i : i in [1..#cs]
+                       | IsConjugate(G,cs[i][3],G!mm)]>;
+     /*                  */
+
+  tl := AssociativeArray();
+
+  for c in cc do
+    mat := phi(c[3]);
+    mm := ExtractBlock(mat,1,1,4,4); 
+    x := <Order(mm), CycleStructure(c[3])@untally>;
+    if not IsDefined(tl, x) then
+      tl[x] := <0,0>;
+    end if;
+    t := ZZ!tr(mat);
+    /* // should tally also using the class in G
+    if x in [[6,6,6,1,1],[6,3,3,3,3,1,1],[4,2,2,2,2,2,2,2,2],[6,6,3,3,2],[6,6,6,2],[10,5,5]] then
+      mm := ExtractBlock(mat,1,1,4,4); 
+      print x, t, Order(mm), #Conjugates(G, G!mm);
+    end if;
+    */
+    /* if tl[x][2-t] gt 0 then */
+    /*   print x; */
+    /* end if; */
+    tl[x][t+1] +:= c[2];
+    /* if tr(phi(c[3])) eq 0 then */
+    /*   tl[x] := <tl[x][1]+c[2], tl[x][2]>; */
+    /* else */
+    /*   tl[x] := <tl[x][1], tl[x][2]+c[2]>; */
+    /* end if; */
+  end for;
+  
+  return tl;
+  //[<x,tl[x]> :  x in Sort(Setseq(Keys(tl)))];
+  /* Sort([ */ 
+  /*   <CycleStructure(c[3])@untally, tr(phi(c[3]))> */
+  /* : c in cc]); */
+end function;
+
+
+
+
+intrinsic ModPrimes(m :: Modularity, R2 :: FldAb, phi :: Map
+                      : debug := false) -> []
+  { .. }
+  F := AbsoluteField(NumberField(R2));
+  //f2 := DefiningPolynomial(F);
+  ZF := MaximalOrder(F);
+  ZR1 := MaximalOrder(m`R1);
+  tl := show_classes(m, phi);
+  //print [<a,tl[a]> : a in Keys(tl)];
+  if &and [tl[a,2] eq 0 : a in Keys(tl)] then
+    return 1;
+  end if;
+  for p in PrimesUpTo(10000) do
+    if p in [2,277] then continue; end if; // CHECK DISC OF R1 !!!
+    //pat := factpat(f2, p);
+    pat10 := degrees(Factorization(p*ZR1));
+    ord10 := LCM(pat10);
+    pat := <ord10,degrees(Factorization(p*ZF))>;
+    if debug then printf "%o, %o, ", p, pat; end if;
+    if tl[pat][1] eq 0 then
+      if debug then print tl[pat], "    ***"; end if;
+      return p;
+    elif tl[pat][2] ne 0 then
+      if debug then print tl[pat], "    ?", pat10; end if;
+    else
+      if debug then print tl[pat]; end if;
+    end if;
+  end for;
+  //return f2;
+  //return [];
+end intrinsic;
+
+

ex-277.m

+// Galois theory for C277
+
+ZZ := IntegerRing();
+ZZx<x> := PolynomialRing(ZZ);
+
+//f10 := x^10 + 6*x^9 + 15*x^8 + 22*x^7 + 24*x^6 + 12*x^5 - 18*x^4 - 34*x^3 - 9*x^2 + 2*x - 1;
+
+f10 := x^10 + 3*x^9 + x^8 - 10*x^7 - 17*x^6 - 7*x^5 + 11*x^4 + 18*x^3 + 13*x^2 + 5*x + 1;
+
+R1 := NumberField(f10);
+
+P2 := Radical(2*Integers(R1));
+S0 := Radical(277*Integers(R1));
+S11 := P2^11 * S0;
+
+gal, r, s := GaloisGroup(f10);
+
+// all the quadratic extensions ramified at (S, inf)
+function quad_ext(S, inf)
+  ray, m := RayClassGroup(S, inf);
+  return
+  [ AbelianExtension(Inverse(mq)*m)
+    where _, mq is quo<ray | Q>
+    where Q is QQ`subgroup
+  : QQ in Subgroups(ray : Quot:=[2])];
+end function;
+
+qexts := quad_ext(S11, [1..#RealPlaces(R1)]);
+
+//load "modularity.m";
+//load "mod-S5b.m";
+
+// m : gal -> G
+_, rho := IsIsomorphic(gal, G1);
+assert rho eq rho;
+
+
+

ex.m

+//AttachSpec("mod/mod.spec");
+
+//////////////////////////////////////////////////
+
+GSp4 := Sp(4,2);
+
+ex_skip := false;
+
+// S5(B)
+if not ex_skip or not assigned(G1) then
+// printf ".";
+G1 := [H`subgroup
+      : H in Subgroups(GSp4 : OrderEqual:=120)
+      | &and [Trace(x) eq 1 : x in order36]
+        where order36 := [x : x in H`subgroup | Order(x) in [3,6]]
+      ];
+assert #G1 eq 1; G1 := G1[1];
+H1 := [H`subgroup
+      : H in Subgroups(G1 : OrderEqual:=12)
+      | #ConjugacyClasses(H`subgroup) eq 6 ];
+assert #H1 eq 1; H1 := H1[1];
+assert IsIsomorphic(H1, DihedralGroup(6));
+assert IsTrivial(&meet Conjugates(G1,H1));
+end if;
+
+// S5(A)
+if not ex_skip or not assigned(G2) then
+// printf ".";
+G2 := [H`subgroup
+      : H in Subgroups(GSp4 : OrderEqual:=120)
+      | &and [Trace(x) eq 0 : x in order36]
+        where order36 := [x : x in H`subgroup | Order(x) in [3,6]]
+      ];
+assert #G2 eq 1; G2 := G2[1];
+H2 := [H`subgroup
+      : H in Subgroups(G2 : OrderEqual:=6)
+      | &+ [#Conjugates(G2, x) : x in H`subgroup] eq 86];
+assert #H2 eq 1; H2 := H2[1];
+assert IsIsomorphic(H2, SymmetricGroup(3));
+assert #Conjugates(G2,ConjugacyClasses(H2)[2][3]) eq 15;
+assert IsTrivial(&meet Conjugates(G2,H2));
+end if;
+
+// S6
+if not ex_skip or not assigned(G3) then
+// printf ".";
+G3 := GSp4;
+H3 := [H`subgroup
+      : H in Subgroups(G3 : OrderEqual:=36)
+      | not IsEmpty(order6)
+        and &and [Trace(x) eq 1 : x in order6]
+        where order6 := [x : x in H`subgroup | Order(x) eq 6]
+      ];
+assert #H3 eq 1; H3 := H3[1];
+assert GroupName(H3) eq "S3^2";
+assert IsTrivial(&meet Conjugates(G3,H3));
+end if;
+
+
+// A6
+if not ex_skip or not assigned(G4) then
+// printf ".";
+G4 := Subgroups(GSp4:OrderEqual:=360);
+assert #G4 eq 1; G4 := G4[1]`subgroup;
+H4 := [H`subgroup
+      : H in Subgroups(G4:OrderEqual:=18)];
+assert #H4 eq 1; H4 := H4[1];
+assert GroupName(H4) eq "C3:S3";
+assert IsTrivial(&meet Conjugates(G4,H4));
+end if;
+
+// S3wrC2
+if not ex_skip or not assigned(G5) then
+// printf ".";
+G5 := Subgroups(GSp4:OrderEqual:=72);
+assert #G5 eq 1; G5 := G5[1]`subgroup;
+H5 := [H`subgroup
+      : H in Subgroups(G5:OrderEqual:=4)
+      | #&meet [ Kernel( hom<gl4->gl4 |
+                [h*a-a*h : a in OrderedGenerators(gl4)]> ) 
+               : h in Generators(H`subgroup) ] eq 64 ]
+        where gl4 is RMatrixSpace(GF(2), 4, 4);
+assert #H5 eq 1; H5 := H5[1];
+assert GroupName(H5) eq "C2^2";
+assert IsTrivial(&meet Conjugates(G5,H5));
+end if;
+
+// A5(B)
+if not ex_skip or not assigned(G6) then
+// printf ".";
+G6 := [H`subgroup
+      : H in Subgroups(GSp4:OrderEqual:=60)
+      | IsAbsolutelyIrreducible(H`subgroup)];
+assert #G6 eq 1; G6 := G6[1];
+H6 := Subgroups(G6:OrderEqual:=6);
+assert #H6 eq 1; H6 := H6[1]`subgroup;
+assert IsTrivial(&meet Conjugates(G6,H6));
+end if;
+
+// S3^2(A)
+if not ex_skip or not assigned(G7) then
+// printf ".";
+G7 := [H`subgroup
+      : H in Subgroups(GSp4 : OrderEqual:=36)
+      | not IsEmpty(order6)
+        and &and [Trace(x) eq 0 : x in order6]
+        where order6 := [x : x in H`subgroup | Order(x) eq 6]
+      ];
+assert #G7 eq 1; G7 := G7[1];
+H7 := [H`subgroup
+      : H in Subgroups(G7:OrderEqual:=2)
+      | #Centralizer(G7, H`subgroup) eq 4];
+assert #H7 eq 1; H7 := H7[1];
+assert IsTrivial(&meet Conjugates(G7,H7));
+end if;
+
+// C3:S3.C2
+if not ex_skip or not assigned(G8) then
+// printf ".";
+G8 := [H`subgroup
+      : H in Subgroups(GSp4 : OrderEqual:=36)
+      | IsEmpty(order6)
+        where order6 := [x : x in H`subgroup | Order(x) eq 6]
+      ];
+assert #G8 eq 1; G8 := G8[1];
+H8 := Subgroups(G8 : OrderEqual:=2);
+assert #H8 eq 1; H8 := H8[1]`subgroup;
+assert IsTrivial(&meet Conjugates(G8,H8));
+end if;
+
+// F5
+if not ex_skip or not assigned(G9) then
+// printf ".";
+G9 := Subgroups(GSp4 : OrderEqual:=20);
+assert #G9 eq 1; G9 := G9[1]`subgroup;
+H9 := Subgroups(G9 : OrderEqual:=2);
+assert #H9 eq 1; H9 := H9[1]`subgroup;
+assert IsTrivial(&meet Conjugates(G9,H9));
+end if;
+
+//print "";
+
+Gs := [G1, G2, G3, G4, G5, G6, G7, G8, G9];
+Hs := [H1, H2, H3, H4, H5, H6, H7, H8, H9];

run-277.m

+Attach("Modularity.m");
+Attach("GrpExt.m");
+
+load "ex.m";
+load "ex-277.m";
+
+m := ModNew(R1, rho);
+PS := Multiset([]);
+time for i := 1 to #qexts do
+  R2:=qexts[i];
+  phis := ModPhis(m, R2);
+  dim := Integers() ! Log(2, #Domain(phis[1])/#G1);
+  ps := [ModPrimes(m, R2, phi) : phi in phis];
+  for p in ps do 
+      Include(~PS, p);
+  end for;
+  printf "%4o: dim =%2o, %o\n", i, dim, Sort(ps);
+
+  if i mod 10 eq 0 then
+      printf "\n%o\n\n", PS;
+  end if;
+
+end for;
+
+printf "\n\nFINAL: %o\n\n", PS;