/*
    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)
   
    This module implements functions that can be used to distinguishi locally conjugate mod-m Galois images that are non-conjugate in GL(2,Z/mZ),
    as discussed in section 5.5 of the paper (the paper only considers prime m, but the implementation below handles any integer m>1)
    
    The function CurveTorsionDegree(E,m) takes an elliptic curve E over a number field K and an integer m and returns the degree of the
    minimal extension L/K over which the base change of E has a rational point of order m.  In the case that m is prime and E admits a rational
    isogeny of degree m, it uses Elkies' algorithm to compute the isogeny kernel polynomial, which speeds things up dramatically.
    
    The function GroupTrosionDegree(G) takes a subgroup of GL(2,Z/mZ) and returns the index of the largest subgroup of G that fixes
    a vector of order m.   If m is prime it will do this efficiently, otherwise it uses brute force.
*/

// Given E1 in short Weierstrass form, a prime ell, instantiated modular polynomials Phi_ell(j1,Y), ((d/dx)Phi_ell)(j1,Y), ((d/dx)^2Phi_ell)(j1,Y), and isogenous j2,
// computes normalized E2 and isogeny kernel poly h(x) of degree (ell-1)/2 whose roots are x-coords of non-trivial kernel elements
// Based on Elkies' algorithm (see his paper in the Atkin birthday conference proceedings MR1486831), does not handle the case j1==j2
IsogenyKernelPolynomial := function(E1,ell,phi,phi_x,phi_xx,j2)
    local a1,b1,a2,b2,c,j1,F,R,phix,phiy,phixx,phixy,phiyy,m1,j1p,j2p,n2,k1,k2,z,p1,d,p,s,t,h;
    c:=Coefficients(E1);
    error if c[1] ne 0 or c[2] ne 0 or c[3] ne 0, "ElkiesKernelPolynomial requires E1 to be in short Weierstrass form y^2=x^2+Ax+B";
    error if ell le 2 or not IsPrime(ell), "ElkiesKernelPolynomial requires prime ell";
    a1:=c[4]; b1:=c[5];
    j1:=jInvariant(E1);
    error if j1 eq 0 or j1 eq 1728, "ElkiesKernelPolynomial requires j(E1) != 0,1728";
    error if j2 eq 0 or j2 eq 1728, "ElkiesKernelPolynomial requires j2 != 0,1728";
    error if j1 eq j2, "ElkiesKernelPolynomial requires j1 != j2";
    F:=Parent(j2);
    phix:=F!Evaluate(phi_x,j2);
    phiy:=F!Evaluate(Derivative(phi),j2);
    phixx:=F!Evaluate(phi_xx,j2);
    phiyy:=F!Evaluate(Derivative(Derivative(phi)),j2);
    phixy:=F!Evaluate(Derivative(phi_x),j2);
    m1:=18*b1/a1;
    j1p:=m1*j1;
    j2p:=-j1p*phix/(ell*phiy);
    m2:=j2p/j2;
    k1:=j1p/(1728-j1);
    k2:=j2p/(1728-j2);
    a2:=ell^4/48*m2*k2;
    b2:=ell^6/864*m2^2*k2;
    z:=-(j1p^2*phixx+2*ell*j1p*j2p*phixy+ell^2*j2p^2*phiyy) / (j1p*phix);
    p1:=ell*(z/2+(k1-ell*k2)/4+(ell*m2-m1)/3);
    d:=ExactQuotient(ell-1,2);
    p:=[F!0];
    p[0+1]:= F!d; p[1+1] := p1/2; p[2+1] := ((1-10*d)*a1-a2) / 30;  p[3+1] := ((1-28*d)*b1-42*p[1+1]*a1-b2) / 70;
    c:=[F!0];
    c[0+1]:=F!0; c[1+1]:=6*p[2+1]+2*a1*p[0+1]; c[2+1]:=10*p[3+1]+6*a1*p[1+1]+4*b1*p[0+1];
    for r:=2 to d-1 do
        t:=&+[c[n+1]*c[r-n+1]:n in [1..r-1]];
        c[r+1+1]:= ((3*t - (2*r-1)*(r-1)*a1*c[r-1+1] - (2*r-2)*(r-2)*b1*c[r-2+1] ) / ((r-1)*(2*r+5)));
    end for;
    for n:=3 to d-1 do
        p[n+1+1] := ( c[n+1] - (4*n-2)*a1*p[n-1+1] - (4*n-4)*b1*p[n-2+1] ) / (4*n+2);
    end for;
    s:=[F!1];
    for i:=1 to d do
        t:=&+[(-1)^k*p[k+1]*s[i-k+1]: k in [1..i]];
        s[i+1]:=-t/i;
    end for;
    R<x>:=PolynomialRing(F);
    h:=R!0;
    for i:= 0 to d do h:= h+(-1)^i*s[i+1]*x^(d-i); end for;
    return h, EllipticCurve([a2, b2]);
end function;
    
// Given E in short Weierstrass form, a prime ell, and the j-invariant of an elliptic curve that is ell-isogenous to E,
// Returns the kernel polynomial for the isogeny.
IsogenyKernelPolynomialFromjInvariant := function(E,ell,j2)
    local c,R,R1,R2,Phi,phi,phix,phixx,j1;
    
    c:=Coefficients(E);
    error if c[1] ne 0 or c[2] ne 0 or c[3] ne 0, "IsogenyKernelPolynomialFromjInvariant requires E to be in short Weierstrass form y^2=x^2+Ax+B";
    error if ell le 2 or not IsPrime(ell), "IsogenyKernelPolynomialFromjInvariant requires prime ell";
    R:=BaseRing(E);
    R1<x>:=PolynomialRing(R,1);
    R2<X,Y>:=PolynomialRing(R,2);
    Phi := R2!ClassicalModularPolynomial(ell);
    j1:=jInvariant(E);
    _,phi:=IsUnivariate(Evaluate(Phi,1,j1));
    error if Evaluate(phi,j2) ne 0, "IsogenyKernelPolynomialFromjInvariant requires j2 to be the j-invariant of an ell-isogenous curve";
    _,phix := IsUnivariate(Evaluate(Derivative(Phi,1),1,j1));
    error if Evaluate(phix,j2) eq 0, "IsogenyKernelPolynomialFromjInvariant requires j2 to be the j-invariant of an ell-isogenous curve connected by a unique ell-isogeny";
    _,phixx := IsUnivariate(Evaluate(Derivative(Derivative(Phi,1),1),1,j1));
    return IsogenyKernelPolynomial(E,ell,phi,phix,phixx,j2);
end function;

// Given E in short Weierstrass form and a prime ell, returns a list of kernel polynomials for each rational ell-isogeny admitted by E
// Does not try to handle the case where E is ell-isogenous to itself, or has jInvariant 0 or 1728 (returns an empty list in this case)
IsogenyKernelPolynomials := function(E,ell)
    local c,R,R1,R2,Phi,j,phi,phix,phixx,rs;
    
    c:=Coefficients(E);
    error if c[1] ne 0 or c[2] ne 0 or c[3] ne 0, "IsogenyKernelPolynomials requires E to be in short Weierstrass form y^2=x^2+Ax+B";
    error if ell le 2 or not IsPrime(ell), "IsogenyKernelPolynomials requires prime ell";
    R:=BaseRing(E);
    R1<x>:=PolynomialRing(R,1);
    R2<X,Y>:=PolynomialRing(R,2);
    Phi := R2!ClassicalModularPolynomial(ell);
    j:=jInvariant(E);
    _,phi:=IsUnivariate(Evaluate(Phi,1,j));
    if Evaluate(phi,j) eq 0 then return []; end if;
    rs:=[r[1]: r in Roots(phi) | r[1] ne 0 and r[1] ne 1728 and r[1] ne j];
    if #rs eq 0  then return []; end if;
    _,phix := IsUnivariate(Evaluate(Derivative(Phi,1),1,j));
    _,phixx := IsUnivariate(Evaluate(Derivative(Derivative(Phi,1),1),1,j));
    return [IsogenyKernelPolynomial(E,ell,phi,phix,phixx,r) : r in rs];
end function;

// Checks whether E1 is ell-isogenous to E2
VerifyIsogeny := function(E1,E2,ell)
    local pi,hs;
    assert IsIsomorphic(BaseRing(E2),BaseRing(E1));
    pi:=hom<BaseRing(E2)->BaseRing(E1)|BaseRing(E1).1>;
    E1:=WeierstrassModel(E1);
    hs:=IsogenyKernelPolynomials(E1,ell);
    if #hs eq 0 then return false; end if;
    return &or[IsIsomorphic(IsogenyFromKernel(E1,h),EllipticCurve([pi(c):c in Coefficients(E2)])):h in hs];
end function;

// Returns a list of curves ell-isogenous to E1
IsogenousCurves := function(E1,ell)
    E1:=WeierstrassModel(E1);
    return [IsogenyFromKernel(E1,h): h in IsogenyKernelPolynomials(E1,ell)];
end function;

// Given E/k returns polynomial in k[x] whose roots are x-coords of points of order m by removing factors of the division polynomial corresponding to divisors of m
PrimitiveDivisionPolynomial := function(E,m)
    local f;
    f:=DivisionPolynomial(E,m);
    for d in Divisors(m) do if d gt 1 and d lt m then f := ExactQuotient(f,$$(E,d)); end if; end for;
    return f;
end function;

// Returns the degree of the minimal extension of the  base field of E over which E acquires a rational point of order m
// When m is prime and E admits a rational (not necessarily unique) isogeny of degree m, it uses the associated kernel polynomial to speed things up.
// The optimization to determine the degree of the splitting field of the kernel polynomial by reducing modulo random primes is not implemented here (factoring the kernel polynomial is already pretty quick)
CurveTorsionDegree:=function(E,m)
    local f,h,R,A,a,d;
    E:=WeierstrassModel(E);
    R<x>:=PolynomialRing(BaseRing(E));
    f,h:=HyperellipticPolynomials(E);
    assert h eq 0;
    // if m is prime and curve admits an m-isogeny, use the kernel polynomial to speed things up
    if m gt 2 and IsPrime(m) and m lt 60 then
        H:=IsogenyKernelPolynomials(E,m);
        if #H gt 0 then
            d:=(m-1);
            for h in H do
                for a in Factorization(h) do
                    if Degree(a[1]) lt d then
                        RR:=quo<R|a[1]>;
                        if IsSquare(RR!f) then d:= Degree(a[1]);
                        elif 2*Degree(a[1]) lt d then d:=2*Degree(a[1]); end if;
                    end if;
                end for;
            end for;
            return d;
        end if;
    end if;
    phi:=PrimitiveDivisionPolynomial(E,m);
    A:=Factorization(phi);
    d:=2*Degree(A[1][1]);
    for a in A do
        if Degree(a[1]) lt d then
            RR:=quo<R|a[1]>;
            if IsSquare(RR!f) then d:= Degree(a[1]);
            elif 2*Degree(a[1]) lt d then d:=2*Degree(a[1]); end if;
        end if;
    end for;
    return d;
end function;

// return a maximal set of pairwise independent column vectors of order m
// (length is psi(m)=m*prod_{p|m}(1+1/p)
CyclicSubspaces:=function(m)
    local M;
    M:=RSpace(Integers(m),2);
    return {sub<M|[i,j]>:i in Divisors(m),j in [0..m-1]|GCD(i,j) eq 1};
end function;

// Compute the minimal index of a subgroup H of GL(2,Z/mZ) that fixes a vector of order m
// Note that Magma wants matrices to act on row vectors from the right, so we transpose the input
// to remain consistent with our convention that matrices act on column vectors from the left.
GroupTorsionDegree:=function(H)
    local m;
    if #Generators(H) eq 0 then return 1; end if;   // handle trivial case separately
    m:=#BaseRing(H);
    H:=sub<GL(2,Integers(m))|[Transpose(h):h in Generators(H)]>;
    return Min([#Orbit(H,Generators(v)):v in CyclicSubspaces(m)]);
end function;

// compute the degree of the minimal extension over which E admits a rational cyclic m-isogeny
CurveIsogenyDegree:=function(E,m)
    local R,b,phi;
    assert jInvariant(E) ne 0 and jInvariant(E) ne 1728;
    R:=PolynomialRing(BaseRing(E),2);
    b,phi:=IsUnivariate(Evaluate(R!ClassicalModularPolynomial(m),1,jInvariant(E)));
    return Degree(Factorization(phi)[1][1]);
end function;

// determine the index of the largest subgroup of H that fixes a cyclic order m submodule of  Z/mZ x Z/mZ
// Note that Magma wants matrices to act on row vectors from the right, so we transpose the input
// to remain consistent with our convention that matrices act on column vectors from the left.
GroupIsogenyDegree:=function(H)
    local m;
    m:=#BaseRing(H);
    H:=sub<GL(2,Integers(m))|[Transpose(h):h in Generators(H)]>;
    return ExactQuotient(#H,Max([#Stabilizer(H,v):v in CyclicSubspaces(m)]));
end function;
