/*
    Magma code related to the paper :
        "Computing images of Galois representations attached to elliptic curves", Forum of Mathematics, Sigma 4 (2016) e4 (79 pages), http://dx.doi.org/10.1017/fms.2015.33
    Copyright (c) 2015 by Andrew V. Sutherland
    You are welcome to use and/or modify this code under GPL version 2 or any later version (but please be sure to cite the paper if you use this code in your research!)
*/

load "subgroups.m";
load "distinguish.m";

det:=function(G) return {Determinant(g):g in G}; end function;
chi:=function(g) return KroneckerSymbol(Integers()!(Trace(g)^2-4*Determinant(g)),Characteristic(BaseRing(g))); end function;
rk:=function(g) M:=MatrixRing(BaseRing(g),2); return 2-Rank(M!g-M![1,0,0,1]); end function;
sig:=function(g) return <Determinant(g),Trace(g),rk(g)>; end function;
Sig:=function(G) return {sig(g):g in G}; end function;
tzr:=function(G) return #{g:g in G|Trace(g) eq 0}/#G; end function;

tchi:=function(t) return KroneckerSymbol(Integers()!(t[2]^2-4*t[1]),Characteristic(Parent(t[1]))); end function;
tpord:=function(t)
    o:=1; a0:=2; a1:=t[2]; d:=t[1]; t1:=d; t2:=a1;
    while 4*d ne a1^2 do a2:=t2*a1-a0*t1; a0:=a1; a1:=a2; d:=t1*d; o+:=1; end while;
    return o,Order(a1/2);
end function;
tpord1:=function(t)
    a0:=2; a1:=t[2]; d:=t[1]; t1:=d; t2:=a1;
    while d ne 1 do a2:=t2*a1-a0*t1; a0:=a1; a1:=a2; d:=t1*d; end while;
    o:=1; t2:=a1; a0:=2;
    while a1^2 ne 4 do a2:=t2*a1-a0; a0:=a1; a1:=a2; o+:=1; end while;
    return o,Order(a1/2);
end function;

siggen:=function(t)
    F:=Parent(t[1]);
    if t eq <F!1,F!2,1> then return {<F!1,F!2,1>,<F!1,F!2,2>}; end if;
    S:={<F!1,F!2,2>};
    a0:=2; a1:=t[2]; d:=t[1]; t1:=d; t2:=a1;
    while d ne 1 or a1 ne 2 do
        Include(~S,<d,a1,d+1 eq a1 select 1 else 0>);
        a2:=t2*a1-a0*t1; a0:=a1; a1:=a2; d:=t1*d;
    end while;
    return S;
end function;

tord:=function(t) m,n:=tpord(t); return m*n; end function; // ignores ell-part of the order
tmat:=function(t)
    p:=Characteristic(Parent(t[1]));
    if tchi(t) eq 0 then return GL(2,p)![t[2]/2,0,0,t[2]/2]; end if;
    if tchi(t) eq 1 then
        R<x>:=PolynomialRing(GF(p));
        r:=[a[1]:a in Roots(x^2-t[2]*x+t[1])];
        return GL(2,p)![r[1],0,0,r[2]];
    end if;
    x:=t[2]/2;
    e:=GF(p)!PrimitiveRoot(p);
    y:=Sqrt((x^2-t[1])/e);
    return GL(2,p)![x,e*y,y,x];
end function;
    
IsLocallyConjugate := function(G,H1,H2)
    p:=Characteristic(BaseRing(G));
    if #H1 ne #H2 then return false; end if;
    if IsConjugate(G,H1,H2) then return true; end if;
    if #H1 mod p ne 0 or #H2 mod p ne 0 then return false; end if;
    if Min({chi(g):g in H1}) lt 0 or Min({chi(g):g in H2}) lt 0 then return false; end if;
    if Sig(H1) ne Sig(H2) or tzr(H1) ne tzr(H2) then return false; end if;
    b,pi:=IsIsomorphic(H1,H2);
    if not b then return false; end if;
    r:=PrimitiveRoot(p);
    B:=sub<GL(2,p)|[1,1,0,1],[1,0,0,r],[r,0,0,1]>;
    S1:=[x:x in Conjugates(G,H1)|x subset B];
    if #S1 eq 0 then return false; end if;
    G1 := S1[1];
    S2:=[x:x in Conjugates(G,H2)|x subset B];
    if #S2 eq 0 then return false; end if;
    G2 := S2[1];
    G1:={[g[1][1],g[2][2]]:g in G1|g[1][2] eq 0};
    G2:={[g[1][1],g[2][2]]:g in G2|g[1][2] eq 0};
    G3:={[g[2],g[1]]:g in G2};
    if G1 eq G3 then return true; end if;
    return false;
end function;
    
// Algorithm 1: compute generators for subgroup with a given signature and trace zero ratio
gfs :=function(s,zr)
    p :=Characteristic(Parent(Random(s)[1]));
    F:=GF(p);
    R<x>:=PolynomialRing(F);
    z:=F!PrimitiveRoot(p);
    GG := GL(2,p);
    // ell = 2
    if p eq 2 then
        if #s eq 1 then return sub<GG|Identity(GG)>; end if;
        if #s eq 2 and <1,0,1> in s then return sub<GG|[0,1,1,0]>; end if;
        if #s eq 2 then return sub<GG|[1,1,1,0]>; end if;
        return GG;
    end if;
    // order divisible by ell cases
    if <1,2,1> in s then
        a := [n:n in GF(p)|not IsSquare(-n)][1];
        d := z^ExactQuotient(p-1,#{t[1]:t in s});
        if #{t:t in s|tchi(t) eq -1} gt 0 then
            return sub<GG|[1,1,0,1],[0,-a,1/a,0],[1,0,0,d]>;
        else
            o:=Max({Order(Roots(x^2-t[2]*x+1)[1][1]):t in s | t[1] eq 1});
            r:= z^ExactQuotient(p-1,o);
            o:=Max({Order(t[1]):t in s});
            d:=z^ExactQuotient(p-1,o);
            t:=[t:t in s|t[1] eq d][1];
            a:=[b[1]: b in Roots(x^2-t[2]*x+t[1])];
            if #a eq 1 then a:=Append(a,a[1]); end if;
            return sub<GG|[1,1,0,1],[r,0,0,1/r],[a[1],0,0,a[2]]>;
        end if;
    end if;
    // Exceptional cases -- note that we know the order is not divisible by ell at this point
    S:= {tpord(t):t in s};
    dets := {t[1]:t in s};
    Z:=sub<GG|[tmat(t):t in s|tchi(t) eq 0]>;
    // A4
    if S eq {1,2,3} and zr lt 3/8 then
        g:=tmat([t:t in s|tpord(t) eq 3][1]);
        h:=tmat([t:t in s|tpord(t) eq 2][1]);
        hs:=[x:x in Conjugates(GG,h)|sig(x*g) in s and sig(x^-1*g) in s];
        for h in hs do H:=sub<GG|g,h,Z>; if Sig(H) eq s and tzr(H) eq 1/4 then return H; end if; end for;
        assert false;
    end if;
    // S4
    if S eq {1,2,3,4} then
        g:=tmat([t:t in s|tpord(t) eq 3][1]);
        h:=tmat([t:t in s|tpord(t) eq 4][1]);
        hs:=[x:x in Conjugates(GG,h)|sig(x*g) in s and sig(x^-1*g) in s];
        for h in hs do H:=sub<GG|g,h,Z>; if Sig(H) eq s then return H; end if; end for;
        assert false;
    end if;
    // A5
    if S eq {1,2,3,5} then
        g:=tmat([t:t in s|tpord(t) eq 3][1]);
        h:=tmat([t:t in s|tpord(t) eq 5][1]);
        hs:=[x:x in Conjugates(GG,h)|sig(x*g) in s and sig(x^-1*g) in s];
        for h in hs do H:=sub<GG|g,h,Z>; if Sig(H) eq s then return H; end if; end for;
        assert false;
    end if;
    // trivial
    if #{tchi(t):t in s} eq 1 then return Z; end if;
    m:=Max({tpord(t):t in s});
    g:=[tmat(t):t in s|tpord(t) eq m][1];
    m:=Max({tord(t):t in s|tchi(t) eq 0});
    z:=[tmat(t):t in s|tord(t) eq m and tchi(t) eq 0][1];
    H:=sub<GG|g,z>;
    // cyclic
    if Sig(H) eq s and Abs(tzr(H)-zr) lt 1/8 then return H; end if;
    // dihedral
    e := ExactQuotient(p-1,#{h:h in H|chi(h) eq 0});
    if chi(g) eq 1 then
        gamma := GG![0,1,1,0];
        r:=GG![1,0,0,PrimitiveRoot(p)];
    else
        gamma := GG![1,0,0,-1];
        a:=GF(p^2).1;
        r:=tmat([Norm(a),Trace(a),0]);
    end if;
    G:=sub<GG|H,gamma>;
    if Sig(G) eq s then return G; end if;
    G:= sub<GG|H,gamma*r^e>;
    assert Sig(G) eq s;
    return G;
end function;
    
verify_gfs := function(p)
    GG := GL(2,p);
    time S:=[G`subgroup:G in Subgroups(GG)];
    t0:=Cputime();
    cnt:=0;
    for i in [1..#S] do
        t1:=Cputime();
        G:=S[i];
        H:=gfs(Sig(G),tzr(G));
        //assert Sig(H) eq Sig(G);
        assert IsLocallyConjugate(GG,G,H);
        cnt+:=1;
        t2:=Cputime();
        if t2-t1 gt 1 then print cnt, "of", #S, "for p =", p, "in", t2-t1,"secs"; end if;
    end for;
    print "Verified", cnt, "of", #S, "for p =", p, "in", t2-t0,"sec";
    return true;
end function;

// Step 2 of Algorithm 4
TorsionRank:=function(E,N,ell)
    e := Valuation(N,ell);
    if e lt 2 then return e; end if;
    if #BaseRing(E) mod ell ne 1 then return 1; end if;
    M:=ExactQuotient(N,ell^e);
    O := E!0;
    // Get a point P of order ell
    P:=O;
    while P eq O do P:= M*Random(E); end while;
    a:=1; X:=ell*P; while X ne O do a+:=1; P:= X; X*:=ell; end while;
    if a eq e then return 1; end if;
    // generate random points Q of ell-power order until we either find a basis for the ell-torsion or prove the ell-power torsion is cyclic.
    while true do
        Q:=M*Random(E);
        if Q eq E!0 then continue; end if;
        a:=1; X:=ell*Q; while X ne O do a+:=1; Q:= X; X*:=ell; end while;
        if a eq e then return 1; end if;
        if WeilPairing(P,Q,ell) ne 1 then return 2; end if;
    end while;
end function;

// Algorithm 4
//   E is an elliptic curve over a finite field
//   S is a prime ell or list of primes ell
//   t is optional, it allows the caller to pass in the trace if known
EpSigs:=function(E,S:t:=TraceOfFrobenius(E))
    if Type(S) ne SeqEnum then S:=[S]; end if;
    q:=#BaseRing(E);
    T:=[*<>:ell in S*];
    for i:=1 to #S do
        ell:=S[i];
        if IsDivisibleBy(q,ell) then continue; end if;
        r:=TorsionRank(E,q+1-t,ell);
        T[i]:=<GF(ell)!q, GF(ell)!t, r>;
    end for;
    return T;
end function;
    
// Variant of Algorithm 6 using small primes (rather than random ones in an interval)
// For CM curves we consider the normalizer of a Cartan (or the Cartan itself if the CM is defined over K) to be surjective
PossiblyNonsurjectivePrimes:=function(E,S,B:K:=BaseRing(E))
    D:=Integers()!Norm(Discriminant(E));
    cm,cmD:=HasComplexMultiplication(E);
    fK:=DefiningPolynomial(K);
    d:=Integers()!Discriminant(fK);
    T:=AssociativeArray(S);
    if cm then
        for l in S do T[l] := <false,false,1,1,1,1>; end for;
    else
        for l in S do T[l] := l le 7 select <false,false> else <false,false,false>; end for;
    end if;
    P:=PrimesInInterval(1,B);
    K0:=BaseRing(E);
    OK:=RingOfIntegers(K);
    if K0 eq Rationals() then
        // For curves defined over Q this is much faster than computing traces individually
        Qtraces:=TracesOfFrobenius(E,B);
    else
        // we only optimize base changes for E/Q
        E:=ChangeRing(E,K);
        K0:=K;
    end if;
    if cm then cmcm:= IsSquare(OK!cmD); end if;
    for i:=1 to #P do
        p:=P[i];
        if IsDivisibleBy(D,p) or IsDivisibleBy(d,p) then continue; end if;
        if K0 eq Rationals() then
            // For base changes we only look at degree-1 primes
            if #Roots(ChangeRing(fK,GF(p))) eq 0 then continue; end if;
            ss:=[EpSigs(ChangeRing(E,GF(p)),S:t:=Qtraces[i])];
        else
            ss:=[EpSigs(Reduction(E,pp[1]),S) : pp in Factorization(ideal<OK|p>) | pp[2] eq 1 and Norm(pp[1]) le B];
        end if;
        U:=S;
        for s in ss do
            if #U eq 0 then break; end if;
            for j:=1 to #S do
                if U[j] eq 0 then continue; end if;
                l:=S[j];  sl:=s[j]; r:=T[l];
                if l eq p then continue; end if;
                if cm and IsDivisibleBy(cmD,l) then continue; end if;
                if l eq 2 then
                    // We don't bother thinking about the CM cases for l eq 2, they are handled in GaloisImages
                else
                    t:=tchi(sl);
                    if l le 7 then
                        if t eq -1 then r[1] := true; end if;
                        if cm then
                            if t eq 1 then r[2] := true; end if;
                            m,n:=tpord1(sl);
                            if sl[2] ne 0 then r[3]:=LCM(r[3],m); end if;
                            r[4]:=LCM(r[4],n);
                            m,n:=tpord(sl);
                            if sl[2] ne 0 then r[5]:=LCM(r[5],m); end if;
                            r[6]:=LCM(r[6],n);
                        else
                            if t eq 0 and sl[3] eq 1 then r[2]:=true; end if;
                        end if;
                    elif sl[2] ne 0 then
                        if t eq 1 then r[1] := true; elif t eq -1 then r[2]:= true; end if;
                        if cm then
                            m,n:=tpord1(sl);
                            r[3]:=LCM(r[3],m);
                            r[4]:=LCM(r[4],n);
                            m,n:=tpord(sl);
                            r[5]:=LCM(r[5],m);
                            r[6]:=LCM(r[6],n);
                        else
                            u:=sl[2]^2/sl[1];
                            if not u in [1,2,4] and u^2-3*u+1 ne 0 then r[3]:=true; end if;
                        end if;
                    elif cm then
                        m,n:=tpord1(sl);
                        r[4]:=LCM(r[4],n);
                        m,n:=tpord(sl);
                        r[6]:=LCM(r[6],n);
                        if t eq 1 then r[1] := true; elif t eq -1 then r[2]:= true; end if;
                    end if;
                end if;
                T[l]:=r;
                if l eq 2 then
                    if r[1] and r[2] then U[j]:=0; end if;
                else
                    if cm then
                        if cmcm then
                            if (2*r[3] ge l-1 and r[4] eq 2) or (r[5] ge l-1 and r[6] eq l-1) then U[j]:=0; end if;
                        else
                            if (r[1] and r[2] and 2*r[3] ge l-1 and r[4] eq 2) or (r[1] and r[2] and r[5] ge l-1 and r[6] eq l-1) then U[j]:=0; end if;
                        end if;
                    else
                        if l le 7 then
                            if r[1] and r[2] then U[j]:=0; end if;
                        else
                            if r[1] and r[2] and r[3] then U[j]:=0; end if;
                        end if;
                    end if;
                end if;
            end for;
        end for;
        S:=[l:l in U|l ne 0];
    end for;
    return S;
end function;


// Variant of Algorithm 5
//   E is an elliptic curve over a number field
//   S is a prime ell or list of primes ell
//   B is a bound on the norms of the primes p to check (it may be increased if it appears to be too low)
//   K is an optional field of definition, if different from the field over whih E is defined
//      if K is specified E will be base-changed to K, but the algorithm can take advantage of the fact that E is
//      defined over a proper subfield (currently only relevant when E is defined over Q)
GaloisImages:=function(E,S,B:K:=BaseRing(E))
    if Type(S) ne SeqEnum then S:=[S]; end if;
    if #S eq 0 then return [**]; end if;
    Sort(~S);
    for ell in S do assert IsPrime(ell); end for;
    Sodd := 2 in S select Remove(S,1) else S;
    OK:=K eq Rationals() select Integers() else RingOfIntegers(K);
    if #Sodd ne 0 then
        D:=Integers()!Norm(Discriminant(E));
        fK:=DefiningPolynomial(K);
        d:=Integers()!Discriminant(fK);
        T:=[*{}:ell in Sodd*];
        cnt:=[0:ell in Sodd];
        zcnt:=[0:ell in Sodd];
        P:=PrimesInInterval(1,B);
        K0:=BaseRing(E);
        if K0 eq Rationals() then
            // For curves defined over Q this is much faster than computing traces individually
            Qtraces:=TracesOfFrobenius(E,B);
        else
            K0:=K;
        end if;
        E:=ChangeRing(E,K);
        for i:=1 to #P do
            p:=P[i];
            if IsDivisibleBy(D,p) or IsDivisibleBy(d,p) then continue; end if;
            if K0 eq Rationals() then
                // For base changes we only look at degree-1 primes
                n:=#Roots(ChangeRing(fK,GF(p)));
                if n eq 0 then continue; end if;
                t:=Qtraces[i];
                Ep:=ChangeRing(E,GF(p));
                s:=EpSigs(Ep,Sodd:t:=t);
                for j:=1 to #Sodd do
                    if Sodd[j] eq p then continue; end if;
                    if not s[j] in T[j] then
                        T[j] := T[j] join siggen(s[j]);
                    end if;
                    cnt[j]+:=n;
                    if s[j][2] eq 0 then zcnt[j]+:=n; end if;
                end for;
            else
                A:=Factorization(ideal<OK|p>);
                for a in A do
                    pp:=a[1];
                    if Norm(pp) le B then
                        Ep:=Reduction(E,pp);
                        s:=EpSigs(Ep,Sodd);
                        for j:=1 to #Sodd do
                            if Sodd[j] eq p then continue; end if;
                            if not s[j] in T[j] then
                                T[j] := T[j] join siggen(s[j]);
                            end if;
                            cnt[j]+:=1;
                            if s[j][2] eq 0 then zcnt[j]+:=1; end if;
                        end for;
                    end if;
                end for;
            end if;
        end for;
    end if;
    G := [**];
    for i:=1 to #Sodd do
        z :=zcnt[i]/cnt[i];
        try
            H:=gfs(T[i],z);
            if Sig(H) ne T[i] then print "gfs sig mismatch!!"; print Sodd[i], T[i], z, Sig(H); assert false; end if;
            G[i]:=GL2SubgroupLabel(H);
        catch e
            printf "B=%o is too small, doubling it...\n", B;
            return $$(E,S,2*B);
        end try;
        C:=GL2LocalConjugates(G[i]);
        if #C ne 1 then
            assert #C eq 2;
            C:=[c:c in C];
            TD:=[GroupTorsionDegree(GL2SubgroupFromLabel(c)):c in C];
            if TD[1] ne TD[2] then
                td:=CurveTorsionDegree(E,Sodd[i]);
                assert td in TD;
                if td eq TD[1] then G[i]:=C[1]; else G[i] := C[2]; end if;
            else
                G[i]:=Sprintf("%o or %o", C[1], C[2]);
                R<x>:=PolynomialRing(Rationals());
                // Hack to distinguish known ambiguity for j(E) in {-11^2,-11*131^3} corresponding to 11B.10.4 and 11B.10.5 (resp.) over Q
                if Sodd[i] eq 11 and (jInvariant(E) eq -11^2 or (DefiningPolynomial(K) eq x^2-x-8 and MinimalPolynomial(jInvariant(E)) eq x^2 - 9476266420414*x - 29164271838685291043)) then
                    X:=[c:c in C|GL2SubgroupContains("11B.10.4",c)];
                    if #X eq 1 then
                        G[i]:=X[1];
                    end if;
                end if;
                if Sodd[i] eq 11 and (jInvariant(E) eq -11*131^3 or (DefiningPolynomial(K) eq x^2-x-8 and MinimalPolynomial(jInvariant(E)) eq x^2 - 2050534*x - 1672446203)) then
                    X:=[c:c in C|GL2SubgroupContains("11B.10.5",c)];
                    if #X eq 1 then
                        G[i]:=X[1];
                    end if;
                end if;
            end if;
        end if;
    end for;
    // handle ell=2 here
    if 2 in S then G := [*["2Cs","2B","2Cn","","","2G"][#GaloisGroup(DivisionPolynomial(ChangeRing(E,K),2))]*] cat G; end if;
    cm,cmD := HasComplexMultiplication(E);
    if cm then
        cmcm:= IsSquare(OK!cmD);
        for i:=1 to #S do
            ell:=S[i];
            s := GL2SubgroupLabel(GL2SubgroupFromLabel(G[i]) meet SL(2,ell));
            if ell gt 2 then
                if not IsDivisibleBy(cmD,ell) then
                    t := Sprintf("%o%o%o[%o]",ell, cmcm select "C" else "N", KroneckerSymbol(cmD,ell) eq 1 select "s" else "n", ell-1);
                else
                    t := Sprintf("%oB[%o]",ell, ell-1);
                end if;
            else
                if not IsDivisibleBy(cmD,ell) then
                    // Note that for ell=2 we use 2Ns=2B and 2Nn=2G (these are the correct analogs of the normalizers of 2Cs and 2Cn, 2-extensions by a trace zero element)
                    t := cmcm select (KroneckerSymbol(cmD,ell) eq 1 select "2Cs" else "2Cn") else (KroneckerSymbol(cmD,ell) eq 1 select "2B" else "2G");
                else
                    t:= "2B";
                end if;
            end if;
            if s eq t then G[i] := ""; end if;
        end for;
        G := [s:s in G|#s gt 0];
    end if;
    return G;
end function;









    
