## Code to analyze elliptic curves over Q with 3-isogenies
## Based on work of M. Bhargava, Z. Klagsbrun, R. Lemke Oliver, and A. Shnidman
## Code mostly written by Robert Lemke Oliver, robert.lemke_oliver@tufts.edu
## Some functions written by Ari Shnidman.
## Please contact RJLO if you have questions regarding this code.

#####################################
##
##  RELEVANT FUNCTIONS:
##
##  AnalyzeCurve(E) 
##    Summarizes what can be said about 3-Selmer ranks of the elliptic curve E/Q possessing a 3-isogeny
##
##  AnalyzeSha(E)
##    Summarizes what can be said (if anything) about 3-torsion in Sha(E_s) for twists of E/Q
##
##  AverageRank(E,isog=0)
##    Determines the bound on the average 3-Selmer rank of twists of E/Q,
##    using the isogeny isog if E has more than one 3-isogeny.
##    In this event, the function IndAverageRank(E) uses both isogenies to get better bounds.
##
##  RankZero(E,isog=0)
##  RankOne(E,isog=0)
##    Returns lower bounds on the proportion of twists for which 
##    the 3-Selmer group has rank 0 or 1.  If E has two 3-isogenies,
##    the functions IndRankZero and IndRankOne will yield better results.
##
##  GlobalSelmerRatio(E,isog=0)
##    Computes the global Selmer ratio c(phi) using the isogeny isog
##
##  GlobalSelmerMeasure(E,m,isog=0)
##    Determines the measure of the set of twists for which c(phi_s) = 3^m
##    Use this and Theorem 1.5 to obtain rank 0 and 3-Selmer rank 1 twists
##
##  GenPoly(E,isog=0)
##    Produces the generating polynomial for the global Selmer measures
##
##  IndGenPoly(E)
##    For E with two 3-isognies, computes the two-variable generating polynomial for the joint Selmer measures
##
##  IndGlobalSelmerMeasure(E,m,n)
##    Determines the joint Selmer measure of T_m(\phi_1) \cap T_n(\phi_2)
##
##  IndZeroMeasure(E)
##  IndOneMeasure(E)
##    Determines the measures of the s for which one of c(\phi_{1,s}) and c(\phi_{2,s}) is either 0 or \pm 1 for curves with two 3-isognies
##    This yields improved rank bounds 
##
##  ShaProps(E)
##    For a curve E which is isogenous to a curve with two 3-isogenies,
##    returns a tuple [x1, x2, ..., xr] where each xi is the minimum proportion of twists
##    for which Sha[3] has 3-rank at least 2*i
##
#####################################


print "\nLoaded code based on the paper \"3-Isogeny Selmer groups and ranks of abelian varieties in quadratic twist families over a number field\" by Bhargava, Klagsbrun, Lemke Oliver, and Shnidman."
print "\nIf you're unsure where to get started, run the command ThreeIsogenyHelp()."
print ""

def ThreeIsogenyHelp():
    print "This code includes various functions to study elliptic curves E/Q possessing a rational 3-isogeny."
    print "\nBasic functions:"
    print "\n*) AnalyzeCurve(E), which takes such a curve E and summarizes what the methods of BKLOS can say."
    print "\n*) AverageRank(E), which returns a bound on the average 3-Selmer rank of the twists of E."
    print "\n*) RankZero(E), which returns a lower bound on the proportion of twists with 3-Selmer rank 0."
    print "\n*) RankOne(E), which returns a lower bound on the proportion of twists with 3-Selmer rank 0."
    print "\n*) GenPoly(E), which returns a generating Laurent polynomial for the measures of the sets T_m(\phi)."
    print "\nApart from AnalyzeCurve, each function takes an optional argument isog=0 or 1, which specifies which isogeny to use in the event that E has two independent 3-isogenies."
    print "\nIf E does have two 3-isogenies, prepending the string Ind to these functions (e.g. IndAverageRank(E)) obtains improved results by exploiting both isogenies."
    print "\nThere are other functions too, documented in the code."
    print "\nThank you for your interest!\n"
    

def PrimeIsogenyDegrees(E):
    return [phi.degree() for phi in E.isogenies_prime_degree()]

def cp(E,p,isog=0):
    #Returns local Selmer ratio of E at p   

    E = E.minimal_model()
    phi = E.isogenies_prime_degree(3)[isog]
    F = phi.codomain()
    F = F.minimal_model()
    f = phi.kernel_polynomial()
    phi = E.isogeny(f,F,3)
    c = F.tamagawa_number(p)/E.tamagawa_number(p)
    if p == 3:
        c = c*phi.formal().coefficients()[0]
    return c

def GlobalSelmerRatio(E,isog=0):
    #Returns the global Selmer ratio of E/Q
    
    if len(E.division_polynomial(3).roots())==0:
        print "The elliptic curve must have a 3-isogeny!"
        return 0
    
    baseline = mul([ cp(E,p,isog) for p in (3*E.conductor()).prime_divisors()])
    phi_root = E.isogenies_prime_degree(3)[isog].kernel_polynomial().roots()[0][0]
    _.<t>=PolynomialRing(Rationals())
    Delta=E.defining_polynomial()(x=phi_root,y=t,z=1).discriminant()
    if Delta>0:
        return 1/3*baseline
    else:
    	return baseline
    
def F2(E,isog):
    #Returns the local factor at p = 2	
    
    f = 0
    # 2-adic square classes and densities
    S = [[1,2],[3,2],[5,2],[7,2],[2,1],[6,1],[10,1],[14,1]]
    for s in S:
    	m = cp(E.quadratic_twist(s[0]),2,isog).valuation(3) 
        f+= x^m*s[1]
    return f/12

def Fp(E,p,isog=0):
    #Returns the factor at p for p > 2    

    f = 0
    u = least_quadratic_nonresidue(p)

    # p-adic square classes and densities 	
    S = [[1,p],[u,p],[p,1],[p*u,1]]
    for s in S:
    	m = cp(E.quadratic_twist(s[0]),p,isog).valuation(3)
    	f+= x^m*s[1]
    return f/(2*p+2) 

def Find2(E):
    #Returns the local factor at p = 2 for a curve with independent isogenies
    
    var('x,y')
    f = 0
    # 2-adic square classes and densities
    S = [[1,2],[3,2],[5,2],[7,2],[2,1],[6,1],[10,1],[14,1]]
    for s in S:
    	m = cp(E.quadratic_twist(s[0]),2,0).valuation(3) 
    	n = cp(E.quadratic_twist(s[0]),2,1).valuation(3) 
        f+= x^m*y^n*s[1]
    return f/12
    
def Findp(E,p):
    #Returns the factor at p for p > 2 for a curve with independent isogenies

    var('x,y')

    f = 0
    u = least_quadratic_nonresidue(p)

    # p-adic square classes and densities 	
    S = [[1,p],[u,p],[p,1],[p*u,1]]
    for s in S:
    	m = cp(E.quadratic_twist(s[0]),p,0).valuation(3)
    	n = cp(E.quadratic_twist(s[0]),p,1).valuation(3)
    	f+= x^m*y^n*s[1]
    return f/(2*p+2) 

def GenPoly(E,isog=0):
    #Returns the generating function for the global Selmer ratio 
    
    if len(E.division_polynomial(3).roots())==0:
        print "The elliptic curve must have a 3-isogeny!"
        return 0
    Ep=E.isogenies_prime_degree(3)[isog].codomain()
    
    var('x')
    
    #Start off with the factor at infinity
    F=(1+1/x)/2 
    
    #Now the finite primes
    N=3*E.conductor()
    for p in N.prime_divisors():
	if p == 2:
	   F*=F2(E,isog)
	if p > 2:
	   F*=Fp(E,p,isog)
    return F.expand()

def IndGenPoly(E):
    if len(E.division_polynomial(3).roots())<2:
        print "The elliptic curve must have two 3-isogenies!"
        return 0
    
    var('x,y')
    
    #Start off with the factor at infinity
    ker0=E.isogenies_prime_degree(3)[0].kernel_polynomial().roots()[0][0]
    ker1=E.isogenies_prime_degree(3)[1].kernel_polynomial().roots()[0][0]
    F=(1/x+1/y)/2
    
    
    #Now the finite primes
    N=3*E.conductor()
    for p in N.prime_divisors():
	if p == 2:
	   F*=Find2(E)
	if p > 2:
	   F*=Findp(E,p)
    return F.expand()

def GlobalSelmerMeasure(E,m,isog=0):
    F=GenPoly(E,isog)
    for C in F.coefficients():
        if C[1] == m:
            return C[0]
    return 0

def AverageRank(E,isog=0):
    var('x')
    F=GenPoly(E,isog)
    bound=add([ C[0]*(abs(C[1])+3^(-abs(C[1]))) for C in F.coefficients()])
    return bound
    
def IndGlobalSelmerMeasure(E,m,n):
    var('x,y')
    F=IndGenPoly(E)
    for C in F.coefficients():
        if C[1] == m:
            A=C[0]
            for a in A.coefficients():
                if a[1]==n:
                    return a[0]
    return 0

def IndOneMeasure(E):
    var('x,y')
    F=IndGenPoly(E)
    G=F(x=y,y=x)
    prop=add([C[0](y=1) for C in F.coefficients() if abs(C[1])==1])
    prop+=add([C[0](y=1) for C in G.coefficients() if abs(C[1])==1])
    return prop-add(add(IndGlobalSelmerMeasure(E,m,n) for m in [-1,1]) for n in [-1,1])

def IndZeroMeasure(E):
    var('x,y')
    F=IndGenPoly(E)
    G=F(x=y,y=x)
    prop=add([C[0](y=1) for C in F.coefficients() if C[1]==0])
    prop+=add([C[0](y=1) for C in G.coefficients() if C[1]==0])
    return prop-IndGlobalSelmerMeasure(E,0,0)

def IndAverageRank(E):
    F=IndGenPoly(E)
    bound=0
    for C in F.coefficients():
        for D in (C[0]*y/y).coefficients():
            m=min(abs(C[1]),abs(D[1]))
            bound += D[0]*(m+3^(-m))
    return bound


    
def AnalyzeCurve(E):
    print "\nAnalyzing the curve with a-invariants",E.ainvs()

    Degs=PrimeIsogenyDegrees(E)
    Has3Iso = false
    for d in Degs:
        if d==3: Has3Iso=true
    
    if len(Degs) == 0:
        print "\nCurve doesn't have any isogenies over Q, so there is nothing to exploit."
        return
    
    if Has3Iso == false:
        print "\nCurve doesn't have 3-isogeny.  Use other methods."
        return
    
    CurveHasTwoIsos = false
    IsoCurveHasTwo = false    
    if len([d for d in Degs if d ==3]) == 2:
    	CurveHasTwoIsos = true
    	print "\nCurve has two 3-isogenies."
    	
    	## Measures of the sets T_0(\phi_0) and T_1(\phi_0)
    	mu00=GlobalSelmerMeasure(E,0)
    	mu01=GlobalSelmerMeasure(E,1)+GlobalSelmerMeasure(E,-1)
    	rb0=AverageRank(E,0)
    	
    	print "\nFrom first isogeny:"
    	print "Average 3-Selmer rank bound:",rb0,"=",RR(rb0)
	print "3-Selmer rank 0 proportion:", (mu00*1/2),"=",RR(mu00*1/2)
    	print "3-Selmer rank 1 proportion:", (mu01*5/6),"=",RR(mu01*5/6)
    	

	## Measures of the sets T_0(\phi_1) and T_1(\phi_1)
    	mu10=GlobalSelmerMeasure(E,0,1)
    	mu11=GlobalSelmerMeasure(E,1,1)+GlobalSelmerMeasure(E,-1,1)
    	rb1=AverageRank(E,1)

    	print "\nFrom second isogeny:"
    	print "Average 3-Selmer rank bound:",rb1,"=",RR(rb1)
	print "3-Selmer rank 0 proportion:", (mu10*1/2),"=",RR(mu10*1/2)
    	print "3-Selmer rank 1 proportion:", (mu11*5/6),"=",RR(mu11*5/6)
    	
    	mu0 = IndZeroMeasure(E)
    	mu1 = IndOneMeasure(E)
    	rb=IndAverageRank(E)
    	
    	print "\nBy combining isogenies:"
    	print "Average 3-Selmer rank bound:",rb,"=",RR(rb)
    	print "3-Selmer rank 0 proportion:", (mu0*1/2),"=",RR(mu0*1/2)
    	print "3-Selmer rank 1 proportion:", (mu1*5/6),"=",RR(mu1*5/6)
    	
    	sha0 = (mu0 - mu00)*1/2 + (mu1 - mu01)*5/6
    	if sha0 > 0:
    	    print "\nThe codomain of the first isogeny has a positive proportion of twists with Sha[3] non-trivial.  Its a-invariants are",E.isogenies_prime_degree(3)[0].codomain().ainvs()
    	    print "The proportion is at least",sha0,"=",RR(sha0)
    	    print "Run AnalyzeSha on the isogenous curve for more info and possibly better results."    	
    	
    	sha1 = (mu0 - mu10)*1/2 + (mu1 - mu11)*5/6
    	if sha1 > 0:
    	    print "\nThe codomain of the second isogeny has a positive proportion of twists with Sha[3] non-trivial.  Its a-invariants are",E.isogenies_prime_degree(3)[1].codomain().ainvs()
    	    print "The proportion is at least",sha1,"=",RR(sha1)
    	    print "Run AnalyzeSha on the isogenous curve for more info and possibly better results."
    	
    if len(E.isogenies_prime_degree(3)[0].codomain().isogenies_prime_degree(3)) == 2:
        IsoCurveHasTwo = true
        
        print "\nCurve has a 3-isogeny, and the isogenous curve has an extra isogeny to exploit."
        
        ## Measures of the sets T_0(\phi_0) and T_1(\phi_0)
    	mu00=GlobalSelmerMeasure(E,0)
    	mu01=GlobalSelmerMeasure(E,1)+GlobalSelmerMeasure(E,-1)
    	rb0=AverageRank(E)
    	
    	print "\nFor the given curve:"
    	print "Average 3-Selmer rank bound:", rb0,"=",RR(rb0)
	print "3-Selmer rank 0 proportion:", (mu00*1/2),"=",RR(mu00*1/2)
    	print "3-Selmer rank 1 proportion:", (mu01*5/6),"=",RR(mu01*5/6)
    	
    	mu0=IndZeroMeasure(E.isogenies_prime_degree(3)[0].codomain())
    	mu1=IndOneMeasure(E.isogenies_prime_degree(3)[0].codomain())
    	rb=IndAverageRank(E.isogenies_prime_degree(3)[0].codomain())
    	
    	print "\nFor the isogenous curve:"
    	print "Average 3-Selmer rank bound:",rb,"=",RR(rb)
    	print "3-Selmer rank 0 proportion:", (mu0*1/2),"=",RR(mu0*1/2)
    	print "3-Selmer rank 1 proportion:", (mu1*5/6),"=",RR(mu1*5/6)
    	
    	sha0 = (mu0 - mu00)*1/2 + (mu1 - mu01)*5/6
    	if sha0 > 0:
    	    print "\nThe curve has a positive proportion of twists with Sha[3] non-trivial."
    	    print "The proportion is at least",sha0,"=",RR(sha0)
    	    print "Run AnalyzeSha for more info and possibly better results."
    	
    
    if CurveHasTwoIsos == false and IsoCurveHasTwo == false:
	    print "\nCurve has a single 3-isogeny and no other 3-level structure to exploit."
	    print "Results for the isogenous curve will be the same as for this curve."
	    print "\nAverage 3-Selmer rank bound:",AverageRank(E),"=",RR(AverageRank(E))
	    print "3-Selmer rank 0 proportion:", GlobalSelmerMeasure(E,0)*1/2,"=",RR(GlobalSelmerMeasure(E,0)*1/2)
	    print "3-Selmer rank 1 proportion:", (GlobalSelmerMeasure(E,1)+GlobalSelmerMeasure(E,-1))*5/6,"=",RR((GlobalSelmerMeasure(E,1)+GlobalSelmerMeasure(E,-1))*5/6)
    
    print ""
    
    return

def AnalyzeSha(E):

    print "\nLooking for non-trivial Sha[3] in twists of the curve with a-invariants",E.ainvs()

    Degs=PrimeIsogenyDegrees(E)
    Has3Iso = false
    for d in Degs:
        if d==3: Has3Iso=true
    
    if Has3Iso == false:
    	print "Curve must have a 3-isogeny, and the isogenous curve must have an extra 3-isogeny."
    	return
    
    if len(E.isogenies_prime_degree(3)) == 2:
    	print "This curve has two 3-isogenies.  Run AnalyzeSha on one of the isogenous curves instead."
    	return
    
    if len(E.isogenies_prime_degree(3)[0].codomain().isogenies_prime_degree(3))==1:
    	print "The isogenous curve must have two 3-isogenies."
    	return
    
    Ep=E.isogenies_prime_degree(3)[0].codomain()
    
    var('x,y')
    F=IndGenPoly(Ep)
    if Ep.isogenies_prime_degree(3)[0].codomain().is_isomorphic(E)==false:
    	F = F(x=y,y=x)

    G=GenPoly(E)
    M=0
    for C in G.coefficients():
        if C[1]>M:
            M=C[1]
    shaprops=[0 for i in range(0,M+1)]

    for C in F.coefficients():
        if C[1]<-1:
            for D in (C[0]*y/y).coefficients():
                if abs(D[1]) < abs(C[1]):
                    m=Integers()(abs(C[1]))
                    n=Integers()(abs(D[1]))
                    for i in range(1,(m-n+2)/2):
                        if (n%2)==0:
                            shaprops[i] += (1-(n+3^(-n))/(m-2*i+2))*D[0]
                        if (n%2)==1:
                            shaprops[i] += (1-(n+3^(-n)-1)/(m-2*i+1))*D[0]
    
    if shaprops[1]==0:
        print "Unfortunately, we failed to prove a positive proportion of twists have non-trivial Sha[3]."
        return
    
    if shaprops[2]>0:
        print "\nNote: The following proportions are nested, i.e. the set of twists for which the 3-rank of Sha[3] is at least 4 is included in the set of twists for which it's at least 2."
    
    print ""
    
    for i in range(1,floor((M+3)/2)):
        if shaprops[i]>0:
            print "Proportion with 3-rank of Sha[3] at least",2*i,":",shaprops[i],"=",RR(shaprops[i])
    
    print ""
    
    return


def ShaProps(E):

    Degs=PrimeIsogenyDegrees(E)
    Has3Iso = false
    for d in Degs:
        if d==3: Has3Iso=true
    
    if Has3Iso == false:
    	print "Curve must have a 3-isogeny, and the isogenous curve must have an extra 3-isogeny."
    	return
    
    if len(E.isogenies_prime_degree(3)) == 2:
    	print "This curve has two 3-isogenies.  Run ShaProps on one of the isogenous curves instead."
    	return
    
    if len(E.isogenies_prime_degree(3)[0].codomain().isogenies_prime_degree(3))==1:
    	print "The isogenous curve must have two 3-isogenies."
    	return
    
    Ep=E.isogenies_prime_degree(3)[0].codomain()
    
    var('x,y')
    F=IndGenPoly(Ep)
    if Ep.isogenies_prime_degree(3)[0].codomain().is_isomorphic(E)==false:
    	F = F(x=y,y=x)
    
    G=GenPoly(E)
    M=0
    for C in G.coefficients():
        if C[1]>M:
            M=C[1]
    shaprops=[0 for i in range(0,floor((M+3)/2))]
    
    for C in F.coefficients():
        if C[1]<-1:
            for D in (C[0]*y/y).coefficients():
                if abs(D[1]) < abs(C[1]):
                    m=Integers()(abs(C[1]))
                    n=Integers()(abs(D[1]))
                    for i in range(1,(m-n+2)/2):
                        if (n%2)==0:
                            shaprops[i] += (1-(n+3^(-n))/(m-2*i+2))*D[0]
                        if (n%2)==1:
                            shaprops[i] += (1-(n+3^(-n)-1)/(m-2*i+1))*D[0]
    
    if shaprops[1]==0:
        return [0]
    
    return [s for s in shaprops if s>0]

def RankZero(E,isog=0):
    return 1/2*GlobalSelmerMeasure(E,0,isog)

def RankOne(E,isog=0):
    return 5/6*(GlobalSelmerMeasure(E,1,isog)+GlobalSelmerMeasure(E,-1,isog))

def IndRankZero(E):
    return 1/2*IndZeroMeasure(E)
    
def IndRankOne(E):
    return 5/6*IndOneMeasure(E)