///////////////////////////////////////////////////////////////////////////
// "Growth of torsion groups of elliptic curves upon base changes"
// Enrique González-Jiménez & Filip Najman
///////////////////////////////////////////////////////////////////////////

// 13/10/2016 - Magma 2.21

// Magma script related to Theorem 5.7 in the CM-case

load "subgroups.m";

function Deg(G)
// If the group G acts on the left: M*v (M in G)
// the function returns the cardinal of the orbits 
// excluding (0,0)

  p:=Characteristic(BaseRing(G));
  Gt:=sub<GL(2,p)|[Transpose(g) : g in G]>;
  dv:=Sort([#o : o in Orbits(Gt)]);
  return [dv[k] : k in [2..#dv]];

end function;

function GECM_2(jE)
// Given the j-invariant jE of an elliptic curve E/Q with CM, 
// the function computes a list of possible  G_E(2)
// l=2
// Prop 1.15 Zywina:

  G1:=sub<GL(2,2)|[[1,0,0,1]] >; 
  G2:=sub<GL(2,2)|[[1,0,0,1],[1,0,1,1]] >;
  GG:=GL(2,2);
  if jE in {2^4*3^3*5^3,2^3*3^3*11^3,-3^3*5^3,3^3*5^3*17^3,2^6*5^3} then GE:=[G2]; end if;
  if jE in {-2^(15)*3*5^3,-2^(15),-2^(15)*3^3,-2^(18)*3^3*5^3,-2^(15)*3^3*5^3*11^3,-2^(18)*3^3*5^3*23^3*29^3} then GE:=[GG]; end if;
  if jE eq 1728 then GE:=[G1,G2]; end if;
  if jE eq 0 then GE:=[G2,GG]; end if;
  return GE;

end function;

function DegCM_2(D)
// Let E/Q be a CM elliptic curve by and order of discriminant D, 
// then the function computes the possible degrees [Q(P):Q] where P in E and ord(P)=2
// l=2
// Prop 1.15 Zywina:

  G1:=sub<GL(2,2)|[[1,0,0,1]] >; 
  G2:=sub<GL(2,2)|[[1,0,0,1],[1,0,1,1]] >;
  GG:=GL(2,2);
  if D in {-12,-16,-7,-28,-8} then dv:=[Deg(G2)]; end if;
  if D in {-43,-27,-11,-163,-19,-67} then dv:=[Deg(GG)]; end if;
  if D eq -4 then dv:=[Deg(G1),Deg(G2)]; end if;
  if D eq -3 then dv:=[Deg(G2),Deg(GG)]; end if;
  return dv;

end function;


function DegsCM_jne0(l,D)
// Let E/Q be a CM elliptic curve by and order of discriminant D (\ne -3, that is j(E)\ne 0) and l an odd prime,
// then the function computes the possible degrees [Q(P):Q] where P in E and ord(P)=l
// l odd
// Prop 1.14 Zywina:
// j ne 0
  assert IsPrime(l) and D ne -3 and IsOdd(l);
  Fl:=FiniteField(l);
  
  case IsDivisibleBy(D,l):
    when false: 
      case IsSquare(Fl!D):
        when true: 
                  lNs:=GL2SubgroupFromLabel(IntegerToString(l) cat "Ns");
                  return [Deg(lNs)];
        when false:
                  lNn:=GL2SubgroupFromLabel(IntegerToString(l) cat "Nn");
                  return [Deg(lNn)];
     end case;
   when true:  
     Fl0:=[a : a in Fl | not IsZero(a)];
     Fl2:=[a : a in Fl0 | IsSquare(a)];
     G:=sub<GL(2,l)|[[a,b,0,a] : a in Fl0, b in Fl],[[a,b,0,-a] : a in Fl0, b in Fl] >;
     H1:=sub<GL(2,l)|[[a,b,0,a] : a in Fl2, b in Fl],[[a,b,0,-a] : a in Fl2, b in Fl] >;
     H2:=sub<GL(2,l)|[[a,b,0,a] : a in Fl2, b in Fl],[[-a,b,0,a] : a in Fl2, b in Fl] >;
     return [Deg(H1),Deg(H2),Deg(G)];
  end case;
end function;



function DegsCM_jeq0(l)
// Let E/Q be a CM elliptic curve with j(E)=0 and l an odd prime,
// then the function computes the possible degrees [Q(P):Q] where P in E and ord(P)=l
// l odd

// l odd
// Prop 1.16 Zywina:
// j=0
  assert IsPrime(l) and IsOdd(l);
  lNn:=GL2SubgroupFromLabel(IntegerToString(l) cat "Nn");
  lNs:=GL2SubgroupFromLabel(IntegerToString(l) cat "Ns");

  case l mod 9:
    when 1:           
           return [Deg(lNs)];
    when 8:
           return [Deg(lNn)];  
    when 4:
           Fl:=FiniteField(l);
           Fl0:=[a : a in Fl | not IsZero(a)];
           H:=sub<lNs|[[a,0,0,b] : a in Fl0, b in Fl0 | IsPower(a/b,3)],[[0,a,b,0] : a in Fl0, b in Fl0 | IsPower(a/b,3)] >;
           return [Deg(H),Deg(lNs)];
    when 7:
           Fl:=FiniteField(l);
           Fl0:=[a : a in Fl | not IsZero(a)];
           H:=sub<lNs|[[a,0,0,b] : a in Fl0, b in Fl0 | IsPower(a/b,3)],[[0,a,b,0] : a in Fl0, b in Fl0 | IsPower(a/b,3)] >;
           return [Deg(H),Deg(lNs)];
    when 2:
           lCn:=GL2SubgroupFromLabel(IntegerToString(l) cat "Cn");
           Sub:=Subgroups(lCn: IndexEqual:=3);
           GG:=[x`subgroup : x in Sub];
           assert #GG eq 1;
          H:=sub<GL(2,l)|Generators(GG[1]) join {GL(2,l)![1,0,0,-1]}>;
          return [Deg(H),Deg(lNs)];
    when 5:
           lCn:=GL2SubgroupFromLabel(IntegerToString(l) cat "Cn");
           Sub:=Subgroups(lCn: IndexEqual:=3);
           GG:=[x`subgroup : x in Sub];
           assert #GG eq 1;
          H:=sub<GL(2,l)|Generators(GG[1]) join {GL(2,l)![1,0,0,-1]}>;
          return [Deg(H),Deg(lNs)];                    
    when 3:
          F3:=FiniteField(3);
          F30:=[a : a in F3 | not IsZero(a)];
          H11:=sub<GL(2,3)| [[1,0,0,a] : a in F30]>;
          H31:=sub<GL(2,3)| [[1,b,0,a] : a in F30, b in F3]>;
          H32:=sub<GL(2,3)| [[a,b,0,1] : a in F30, b in F3]>;
          G1:=GL2SubgroupFromLabel(IntegerToString(3) cat "Cs");
          G3:=GL2SubgroupFromLabel(IntegerToString(3) cat "B");
          return [Deg(H11),Deg(H31),Deg(H32),Deg(G1),Deg(G3)]; 
   end case;
end function;

function DegsCM(l,D)
// Let E/Q be a CM elliptic curve by and order of discriminant D and l a prime,
// then the function computes the possible degrees [Q(P):Q] where P in E and ord(P)=l
 if l eq 2 then return DegCM_2(D); end if;
 if D eq -3 then
    return DegsCM_jeq0(l);
   else
    return DegsCM_jne0(l,D);
  end if;
end function;

////////////////////////////////////////////////////////////////

for l in [2,3,5,7,11,13,17,37] do
 printf "\n[l=%o]\n",l;
 Ll:={};
  for D in [-3,-12,-27,-4,-16,-7,-28,-8,-11,-19,-43,-67,-163] do 
    for S in DegsCM(l,D) do
      Ll:=Ll join Seqset(S);
    end for;
  end for;
  print Ll;
end for;

/* OUTPUT: For each prime l below the corresponding degrees [Q(P):Q] where P in E, ord(P)=l and E/Q a CM-elliptic curve


[l=2]
{ 1, 2, 3 }

[l=3]
{ 1, 2, 3, 4, 6, 8 }

[l=5]
{ 8, 16, 24 }

[l=7]
{ 3, 6, 12, 21, 24, 36, 42, 48 }

[l=11]
{ 5, 10, 20, 40, 55, 80, 100, 110, 120 }

[l=13]
{ 24, 48, 96, 144, 168 }

[l=17]
{ 32, 256, 288 }

[l=37]
{ 72, 1296, 1368 }

*/