changes for 353

Grp.m

+intrinsic OuterAutomorphismRepresentatives(G :: Grp) -> []
+  { representatives for Out(G)/Inn(G) }
+  aut := AutomorphismGroup(G);
+  // m1 : autfp -> aut
+  autfp, m1 := FPGroup(aut);
+  // m2 : autfp -> outfp
+  outfp, m2 := OuterFPGroup(aut);
+  // m3 : outfp -> out
+  out, m3 := PermutationGroup(outfp);
+
+  return [g @@ m3 @@ m2 @ m1 : g in out];
+end intrinsic;
+
+intrinsic OuterAutomorphismGenerators(G :: Grp) -> []
+  { generators for Out(G)/Inn(G) }
+  aut := AutomorphismGroup(G);
+  // m1 : autfp -> aut
+  autfp, m1 := FPGroup(aut);
+  // m2 : autfp -> outfp
+  outfp, m2 := OuterFPGroup(aut);
+
+  return [g @@ m2 @ m1 : g in Generators(outfp)];
+end intrinsic;
+
+intrinsic IsIsomorphicPerm(G :: GrpPerm, H :: GrpPerm) -> BoolElt, Map
+  { Assuming G, H have the same degree, return an isomorphism which
+    preserves the cycle decomposition, if such isomorphism exists }
+
+    n := Degree(G);
+    if Degree(H) ne n then
+        return false;
+    end if;
+
+    b, x := IsConjugate(Sym(n), G, H);
+
+    if not b then
+        return false;
+    end if;
+
+    gens := OrderedGenerators(G);
+    isom := Isomorphism(G, H, gens, [x^-1*g*x : g in gens]);
+
+    return true, isom;
+end intrinsic;
+

Modularity.m

@@ -3,20 +3,25 @@
 ///////////////////////
 
 Attach("GrpExt.m");
+Attach("Grp.m");
 
 declare type Modularity;
 
 declare attributes Modularity :
-  R1, VGs, VHs, Vs, auts;
+  R1, gal, rho, 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
+intrinsic ModNew (R1 :: FldNum, rho :: Map) -> Modularity
   { precomputation for modularity test }
+  gal, r, s := GaloisGroup(R1);
+  assert (gal eq Domain(rho));
+  assert GaloisProof(R1, s);
+
   G := Image(rho);
-  H := rho(Stabilizer(Domain(rho), 1));
+  H := rho(Stabilizer(gal, 1));
   SG := SmallParabolic(G);
   VGs := FilterConjugates(Subextensions(SG));
   auts := AssociativeArray();
@@ -31,6 +36,8 @@ intrinsic ModNew (R1 :: FldNum, rho :: HomGrp) -> Modularity
 
   m := New(Modularity);
   m`R1 := R1;
+  m`gal := gal;
+  m`rho := rho;
   m`VGs := VGs;
   m`VHs := VHs;
   m`Vs := Vs;
@@ -39,38 +46,213 @@ intrinsic ModNew (R1 :: FldNum, rho :: HomGrp) -> Modularity
   return m;
 end intrinsic;
 
+intrinsic galois_projection(R1 :: FldNum[FldRat], R2 :: FldAb) -> Map
+  { .. }
+  return galois_projection(R1, NumberField(R2));
+end intrinsic;
+
+intrinsic galois_projection(R1 :: FldNum[FldRat], R2 :: FldNum[FldNum]) -> Map
+  { .. }
+  return galois_projection(R1, AbsoluteField(R2));
+end intrinsic;
+
+
+// Assume R2/R1 is quadratic (?)
+// Find a good map from Gal(R2) --> Gal(R1)
+intrinsic galois_projection(R1 :: FldNum[FldRat],
+                            R2 :: FldNum[FldRat]) -> Map
+  { .. }
+  assert R1 subset R2;
+  assert Degree(R2) eq 2 * Degree(R1);
+
+  gal, r, s := GaloisGroup(R1);
+  assert GaloisProof(R1, s);
+
+  gal2, r2, s2 := GaloisGroup(R2);
+  assert GaloisProof(R2, s2);
+
+  // subgroup of gal2 corresponding to R2
+  gal_R2 := Stabilizer(gal2, 1);
+
+  // the fixed points should be two roots in R2
+  // which are conjugates over R1
+  fp := Exclude(Fix(gal_R2), 1);
+
+  for p in fp do
+
+    // subgroup of gal2 corresponding to R1
+    gal_R1 := Stabilizer(gal2, {1,p});
+    assert #gal_R1 eq 2*#gal_R2;
+
+    // subgroup of gal2 corresponding to R (Galois closure of R1)
+    gal_R := &meet Conjugates(gal2, gal_R1);
+
+    // Check that {1, p} are really conjugates
+    if #gal2 eq #gal_R * #gal then
+      break;
+    end if;
+
+  end for;
+
+  // proj : gal2 --> gal1 
+  proj, gal1 := CosetAction(gal2, gal_R1);
+
+  // isom : gal1 --> gal
+  b, isom := IsIsomorphicPerm(gal1, gal);
+
+  assert b;
+  
+  assert &and [ CycleStructure(c[3]) eq CycleStructure(isom(c[3]))
+              : c in ConjugacyClasses(gal1) ];
+
+  return proj * isom;
+
+end intrinsic;
+
+intrinsic find_isom( proj :: Map, rho :: Map, VG :: GrpExt ) -> Map
+  { find one isomorphism between the domains that fixes the image }
+
+  G := Group(VG);
+  assert Image(rho) eq G;
+  assert #Domain(rho) eq #G;
+
+  VG1 := Domain(proj);
+  VG2 := Self(VG);
+
+  b, phi := IsIsomorphic(VG1, VG2);
+
+  if not b then
+    return false;
+  end if;
+
+  V1 := Kernel(proj);
+  V2 := Abelian(VG);
+
+  // Is it always true?
+  assert phi(V1) eq V2;
+
+  aut := AutomorphismGroup(VG`self);
+  // m1 : autfp -> aut
+  // m2 : autfp -> autp
+  autfp, m1 := FPGroup(aut);
+  autp, m2 := PermutationGroup(autfp);
+  auts := [ psi @@ m2 @ m1 : psi in autp ];
+
+  for a in auts do
+    ok := &and [ VG`pi(a(phi(g))) eq rho(proj(g))
+                 : g in Generators(VG1)];
+
+    if ok then
+      return phi*a;
+    end if;
+  end for;
+
+  return false;
+
+end intrinsic;
+
 intrinsic ModPhis(m :: Modularity, R2 :: FldAb) -> []
   { return a list of all possible phi }
+  return ModPhis(m, NumberField(R2));
+end intrinsic;
+
+
+intrinsic ModPhis(m :: Modularity, R2 :: FldNum) -> Map, []
+  { return proj:gal2->gal1 and a list of all possible phi:gal2->VG }
 
-  assert NumberField(BaseRing(R2)) eq m`R1;
+  proj := galois_projection(m`R1, R2);
 
-  // may not be the most efficient way to do this
-  //gal, r, s := AbsoluteGaloisGroup(R2);
-  gal, r, s := GaloisGroup(AbsoluteField(NumberField(R2)));
+  gal2 := Domain(proj);
 
   ans := [* *];
 
   for VG in m`VGs do
-    b, phi := IsIsomorphic(gal, VG`self);
-    VH := m`VHs[VG]`self;
-    V := m`Vs[VG];
+    b, phi := IsIsomorphic(gal2, VG`self);
     if b then
-      //print Dimension(VG);
+      //print "dim:", 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;
+  return proj, ans;
+
+end intrinsic;
+
+intrinsic ModPhis1(m :: Modularity, R2 :: FldNum) -> Map, []
+  { return proj:gal2->gal1 and a list of all possible phi:gal2->VG }
+
+  proj := galois_projection(m`R1, R2);
+
+  gal2 := Domain(proj);
+
+  ans := [* *];
+
+  ans1 := [* *];
+
+  i := 0;
+  for VG in m`VGs do
+    // phi : gal2 --> VG
+    b, phi := IsIsomorphic(gal2, VG`self);
+    if b then
+      i := i + 1;
+      print i, "dim:", Dimension(VG);
+
+      // pi : VG --> G;
+      pi := VG`pi;
+
+      // phi * pi == proy * isom * rho ?
+      assert Domain(phi) eq gal2;
+      assert Codomain(phi) eq Domain(pi);
+      assert Image(pi) eq Image(m`rho);
+
+      VH := m`VHs[VG]`self;
+      V := m`Vs[VG];
+      //print IsIsomorphic(NumberField(GaloisSubgroup(s, VH@@phi)), m`R1);
+      //print Dimension(VG);
+      
+      // autmap : l
+      aut := AutomorphismGroup(VG`self);
+      // m1 : autfp -> aut
+      // m2 : autfp -> autp
+      autfp, m1 := FPGroup(aut);
+      autp, m2 := PermutationGroup(autfp);
+      auts := [ psi @@ m2 @ m1 : psi in autp ];
+
+      rho := m`rho;
+
+      print "#aut:", #auts;
+
+      for psi in auts do //m`auts[VG] do
+        phi1 := phi * psi;
+        ok := &and [ pi(phi1(g)) eq rho(proj(g))
+                   : g in Generators(gal2)];
+        //if phi1 * pi eq proy * isom * rho then
+        if ok then
+          //print "yes";
+          Append(~ans, phi1);
+        else
+          Append(~ans1, phi1);
+        end if;
+        //assert psi(V) eq V;
+        //if psi(VH) eq VH then // IS IT RIGHT TO FILTER LIKE THIS?
+        //else
+        //  Append(~ans1, phi*psi);
+        //end if;
+      end for;
+      print <#ans, #ans1>;
+    end if;
+  end for;
+
+  print <#ans,#ans1>;
+  return proj, ans;
 
 end intrinsic;
 
 /////////////////////////
 // sorted degrees for a factorization (from high to low)
-degrees := func<fact | Reverse(Sort([Degree(f[1]): f in fact]))>;
+degrees := func<fact | Reverse(Sort([AbsoluteDegree(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
@@ -86,81 +268,116 @@ 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)]>;
-     /*                  */
+function show_classes(pi, phi)
+  gal2 := Domain(pi);
+  assert Domain(phi) eq gal2;
 
   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>;
+  for c in ConjugacyClasses(gal2) do
+
+    x := < CycleStructure(pi(c[3])) @ untally,
+           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; */
+
+    t := ZZ!tr(phi(c[3]));
+
+    tl[x][t+1] +:= 1; //c[2];
   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;
 
 
+function frobenius_classes(pi)
+  gal2 := Domain(pi);
+  ans := AssociativeArray();
 
+  for c in ConjugacyClasses(gal2) do
+
+    x := < CycleStructure(pi(c[3])) @ untally,
+           CycleStructure(c[3])     @ untally >;
+
+    if not IsDefined(ans, x) then
+      ans[x] := [ c[3] ];
+    else
+      Append(~ans[x], c[3]);
+    end if;
+  end for;
+  
+  return ans;
+end function;
 
 intrinsic ModPrimes(m :: Modularity, R2 :: FldAb, phi :: Map
                       : debug := false) -> []
   { .. }
-  F := AbsoluteField(NumberField(R2));
+  return ModPrimes(m, NumberField(R2), phi : debug := debug);
+end intrinsic;
+
+intrinsic AbsoluteDegree( p :: RngOrdIdl) -> FldRatElt
+  { Degree for prime ideals which works for relative extensions}
+  return Factorization(AbsoluteNorm(p))[1,2];
+end intrinsic;
+
+//  pi : Gal(R2) --> Gal(R1)
+// phi : Gal(R2) --> VG
+intrinsic ModPrimes(m :: Modularity, R2 :: FldNum,
+                    proj :: Map, phi :: Map
+          : debug := false, keep_going := false, PrimeBound := 500) -> []
+  { .. }
+  //F := AbsoluteField(R2);
   //f2 := DefiningPolynomial(F);
-  ZF := MaximalOrder(F);
   ZR1 := MaximalOrder(m`R1);
-  tl := show_classes(m, phi);
+  ZR2 := MaximalOrder(R2);
+  discR1 := Discriminant(m`R1);
+  //tl := show_classes(proj, phi);
+  fc := frobenius_classes(proj);
   //print [<a,tl[a]> : a in Keys(tl)];
-  if &and [tl[a,2] eq 0 : a in Keys(tl)] then
+  /*
+  if &and [tl[a,1] ne 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))>;
+  */
+  for p in PrimesUpTo(PrimeBound) do
+    if discR1 mod p eq 0 then
+      continue;
+    end if;
+    pat := < degrees(Factorization(p*ZR1)),
+             degrees(Factorization(p*ZR2)) >;
     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;
+    //if not pat in Keys(tl) then
+    if not pat in Keys(fc) then
+      if debug then print "    not in this extension (?)"; end if;
+      if not keep_going then return p; end if;
+      continue;
+    end if;
+    all_classes := fc[pat];
+    bad_classes := [ c : c in all_classes | tr(phi(c)) eq 0 ];
+    if bad_classes eq [] then
+      if debug then 
+        print "    good prime", [<tr(phi(c)),tr1(phi(c))> : c in all_classes];
+      end if;
+      if not keep_going then return p; end if;
     else
-      if debug then print tl[pat]; end if;
+    if debug then 
+      print [<tr(phi(c)),tr1(phi(c))> : c in all_classes];
     end if;
+   // <#bad_classes, #all_classes-#bad_classes>; end if;
+    end if;
+
+    //elif 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], "    ?"; end if;
+    //else
+    //  if debug then print tl[pat], "     "; end if;
+    //end if;
   end for;
+  return 0; // placeholder
   //return f2;
   //return [];
 end intrinsic;

ex-353.m

+// Galois theory for C353
+
+// faster but conditional
+// SetClassGroupBounds("GRH");
+
+load "ex.m";
+//load "polys.m";
+
+ZZ := IntegerRing();
+ZZx<x> := PolynomialRing(ZZ);
+load "all-353.m";
+f18 := DefiningPolynomial(R1);
+
+gal, r, s := GaloisGroup(f18);
+
+b, rho := IsIsomorphic(gal, G5);
+
+assert b;
+
+// we know Frob5 
+Frob5cycle := 
+ Reverse(Sort([<x,Multiplicity(s, x)> : x in Set(s)]))
+ where s is ([ Degree(f[1]) 
+            : f in Factorization(ChangeRing(f18,GF(5))) ]); 
+
+Frob5 := [c[3] : c in ConjugacyClasses(gal)
+               | CycleStructure(c[3]) eq Frob5cycle ];
+
+assert(#Frob5 eq 1);
+
+Frob5 := Frob5[1];
+
+// we expect Tr(Frob5) = 1
+if Trace(rho(Frob5)) ne 1 then
+  // fix it by composing rho with the outer automorphism
+  rho := rho * a
+  where a is autg`FpGroup[2](outg.1 @@ m)
+  where outg, m is OuterFPGroup(autg)
+  where autg is AutomorphismGroup(G5);
+  print ".";
+end if;
+
+// We expect Tr(Frob5) = 1
+// Make sure this is right
+assert Trace(rho(Frob5)) eq 1;
+
+f353_R1 := GaloisSubgroup(s, H5 @@ rho);
+
+//R1 := NumberField(f353_R1);
+
+P2 := Radical(2*Integers(R1));
+S0 := Radical(353*Integers(R1));
+
+S5 := P2^5 * S0;
+
+//gal, r, s := GaloisGroup(f353_R1);
+
+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(S5, [1..#RealPlaces(R1)]);
+

run-277.m

 Attach("Modularity.m");
 Attach("GrpExt.m");
+Attach("Grp.m");
 
 load "ex.m";
 load "ex-277.m";
@@ -7,10 +8,10 @@ load "ex-277.m";
 m := ModNew(R1, rho);
 PS := Multiset([]);
 time for i := 1 to #qexts do
-  R2:=qexts[i];
-  phis := ModPhis(m, R2);
+  R2:=NumberField(qexts[i]);
+  proj, phis := ModPhis(m, R2);
   dim := Integers() ! Log(2, #Domain(phis[1])/#G1);
-  ps := [ModPrimes(m, R2, phi) : phi in phis];
+  ps := [ModPrimes(m, R2, proj, phi) : phi in phis];
   for p in ps do 
       Include(~PS, p);
   end for;

run-353.m

+Attach("Modularity.m");
+Attach("GrpExt.m");
+Attach("Grp.m");
+
+//load "ex.m";
+load "ex-353.m";
+