Ticket #10973: trac_10973.patch

File trac_10973.patch, 20.9 KB (added by gagansekhon, 11 years ago)
  • sage/schemes/elliptic_curves/all.py

    # HG changeset patch
    # User Author Justin C. Walker <justin@mac.com>
    # Date 1301015058 25200
    # Node ID 673c3b123bf2b49af93407ff533e2fbc473bb67f
    # Parent  ec11b62f0270b890da3b9b2021aed9b30ed718ee
    Add Rado/Anderson code to find integral points on elliptic curves over
    number fields.  This is based on code by John Cremona's student Nook
    (Thotsaphon Thongjunthug).  Minor changes made to handle points on
    curves defined over the rational numbers.
    
    diff -r ec11b62f0270 -r 673c3b123bf2 sage/schemes/elliptic_curves/all.py
    a b  
    4040from ell_curve_isogeny import EllipticCurveIsogeny, isogeny_codomain_from_kernel
    4141
    4242from heegner import heegner_points, heegner_point
     43
     44# Temp until we integrate this more fully:
     45from ell_int_points import *
  • new file sage/schemes/elliptic_curves/ell_int_points.py

    diff -r ec11b62f0270 -r 673c3b123bf2 sage/schemes/elliptic_curves/ell_int_points.py
    - +  
     1#        File: intpts.sage
     2#     Created: Thu Jul 01 04:22 PM 2010 C
     3# Last Change: Fri Jul 02 12:51 PM 2010
     4# Original Magma Code: Thotsaphon "Nook" Thongjunthug
     5# Sage version: Radoslav Kirov, Jackie Anderson
     6
     7from copy import copy
     8from sage.functions.all import sqrt
     9from sage.matrix.all import zero_matrix, MatrixSpace
     10from sage.misc.all import verbose, n, prod
     11from sage.rings.all import is_RealField, polygen, ZZ, RealField, ComplexField
     12from sage.rings.all import QQ, integer_ceil, integer_floor, PolynomialRing
     13
     14# This function should be detached and included as part of projective
     15#  points over number field
     16def abs_log_height(X, gcd_one = True, precision = None):
     17    r''' Computes the height of a point in a projective space over field K.
     18    It assumes the coordinates have gcd equal to 1.
     19    If not use gcd_one flag.
     20    '''
     21    assert isinstance(X,list)
     22    K = X[0].parent()
     23    if precision is None:
     24        RR = RealField()
     25        CC = ComplexField()
     26    else:
     27        RR = RealField(precision)
     28        CC = ComplexField(precision)
     29    places = set([])
     30    if K == QQ:
     31        embs = K.embeddings(RR)
     32        Xideal = X
     33    else:
     34        embs = K.places(precision)
     35        # skipping zero as it currently over K breaks Sage
     36        Xideal = [K.ideal(x) for x in X if not x.is_zero()]
     37    for x in Xideal:
     38        places = places.union(x.denominator().prime_factors())
     39        if not gcd_one:
     40            places = places.union(x.numerator().prime_factors())
     41    if K == QQ:
     42        non_arch = sum([max([RR(P)**(-x.valuation(P)) for x in X]).log() for P in places])
     43    else:
     44        non_arch = sum([P.residue_class_degree() *
     45                        P.absolute_ramification_index() *
     46                        max([x.abs_non_arch(P, precision) for x in X]).log() for P in places])
     47    arch = 0
     48    r,s = K.signature()
     49    for i,f in enumerate(embs):
     50        if i<r:
     51            arch += max([f(x).abs() for x in X]).log()
     52        else:
     53            arch += max([f(x).abs()**2 for x in X]).log()
     54    return (arch+non_arch)/K.degree()
     55   
     56def compute_c6(E,emb):
     57    F = emb.codomain()
     58    if is_RealField(F):
     59        x = polygen(ComplexField(F.prec()))
     60    else:
     61        x = polygen(F)
     62    #f = x**3-27*emb(E.c4())*x-54*emb(E.c6())
     63    f = x**3-emb(E.c4()/48)*x-emb(E.c6()/864)
     64    R = f.roots(multiplicities = False)
     65    m = max([r.abs() for r in R])
     66    return 2*m
     67
     68def compute_c8(L):
     69    # Original code transformed the lattice generators.
     70    # Here we assume they are of standard form.
     71    CC=ComplexField()
     72
     73    w1, w2 = L
     74    m = max(CC(w1).abs(), CC(w2).abs(), CC(w1+w2).abs())
     75    return m
     76
     77# There is a Sage trak ticket implementing this. In the future the library function can be used and this one removed
     78def height_pairing_matrix(points, precision = None):
     79    r = len(points)
     80    if precision is None:
     81        RR = RealField()
     82    else:
     83        RR = RealField(precision)
     84    M = MatrixSpace(RR, r)
     85    mat = M()
     86    for j in range(r):
     87        mat[j,j] = points[j].height(precision = precision)
     88    for j in range(r):
     89        for k in range(j+1,r):
     90            mat[j,k] = ((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2
     91            mat[k,j] = mat[j,k]
     92    return mat
     93
     94def c3(L):
     95   import warnings
     96   with warnings.catch_warnings():
     97       warnings.simplefilter("ignore")
     98       return min(map(abs,height_pairing_matrix(L).eigenvalues()))
     99
     100def h_E(E):
     101    K = E.base_field()
     102    j = E.j_invariant()
     103    c4, c6 = E.c_invariants()
     104    g2, g3 = c4/12, c6/216
     105    return max(abs_log_height([K(1), g2, g3]), abs_log_height([K(1), j]))
     106
     107def h_m(E, P, ElogEmbedP, D7):
     108    K = E.base_field()
     109    M1 =  max([P.height(), h_E(E), D7/K.degree()*abs((ElogEmbedP).log())**2])
     110    M2 =  max([P.height(), h_E(E), D7/K.degree()*abs(ElogEmbedP)**2])
     111#    assert M1==M2
     112    return M2
     113   
     114def Extra_h_m(E, Periods, D7):
     115    D = E.base_field().degree()
     116    h = h_E(E)
     117    return map(lambda x: max([0, h, D7/D*abs(x)**2]), Periods)
     118
     119def d8(E, L, Elog, Periods, D7):
     120    RR=RealField()
     121    C = [RR(1).exp() * h_E(E)]
     122    D = E.base_field().degree()
     123    for i in range(len(L)):
     124        C.append(h_m(E, L[i], Elog[i], D7) / D)
     125    C += [t / D for t in Extra_h_m(E, Periods, D7)]
     126    return max(C);
     127
     128def d9(E, L, Elog, Periods, D7):
     129    RR=RealField()
     130    D = E.base_field().degree()
     131    C = []
     132    for i in range(len(L)):
     133        tmp = RR(1.0).exp() * sqrt(D * h_m(E, L[i], Elog[i], D7) / D7)
     134        C.append( tmp / abs(Elog[i]))
     135    Ehm = Extra_h_m(E, Periods, D7)
     136    C += [RR(1).exp() * sqrt(D*Ehm[i]/D7) / abs(Periods[i]) for i in [0,1]]
     137    verbose("d9: C = %s"%(C,))
     138    return min(C);
     139
     140def d10(E, L, Elog, Periods, D7):
     141    RR=RealField()
     142    D = E.base_field().degree()
     143    n = len(L)+2
     144    scalar = 2 * 10**(8 + 7*n) * (2/RR(1).exp())**(2*n**2)
     145    scalar *= (n+1)**(4*n**2 + 10*n) * D**(2*n + 2)
     146    scalar *= (RR(d9(E, L, Elog, Periods, D7)).log())**(-2*n-1)
     147    for i in range(len(L)):
     148        scalar *= h_m(E, L[i], Elog[i], D7)
     149    scalar *= prod(Extra_h_m(E, Periods, D7))   
     150    return scalar
     151
     152def RHS(D, r, C9, C10, D9, D10, h, Q, expTors):
     153    RR=RealField()
     154    bound = (RR(RR(Q*r*expTors).log()).log() + h + RR(D*D9).log())**(r+3)
     155    bound *= D10*(RR(Q*r*expTors).log() + RR(D*D9).log())
     156    bound += RR(C9*expTors).log()
     157    bound /= C10
     158    return bound
     159
     160def InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors):
     161    RR=RealField()
     162    minQ = max(Q0, RR(D8).exp())
     163    Q = minQ + 1
     164    x = integer_ceil(Q.log(10)) # x = log_10(Q)
     165    verbose("x: %s, %s"%(type(Q), type(x)))
     166    exp_OK = 0   # the exponent that satisfies P.I.
     167    exp_fail = 0 # the exponent that fails P.I.
     168    verbose("Looping on %s < RHS(%s, %s, %s, %s %s, %s, %s, %s, %s"%( 10**(2*x), D, r, C9, C10, D9, D10, h, 10**x, expTors))
     169    while 10**(2*x) < RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
     170        exp_OK = x # Principal Inequality satisfied
     171        x *= 2     # double x, and retry
     172    exp_fail = x # x that fails the Principal Inequality
     173   
     174   
     175    # So now x = log_10(Q) must lie between exp_OK and exp_fail
     176    # Refine x further using "binary search"
     177    while True:
     178        x = integer_floor((exp_OK + exp_fail)/2)
     179        if 10**(2*x) >= RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
     180            exp_fail = x
     181        else:
     182            exp_OK = x
     183        if (exp_fail - exp_OK) <= 1:
     184            break
     185    return exp_fail # over-estimate
     186
     187def Faltings_height(E):
     188    RR=RealField()
     189    K = E.base_field()
     190    c = RR(2).log()
     191    if E.b2().is_zero():
     192        c = 0
     193    h1 = abs_log_height([K(E.discriminant()), K(1)])/6
     194    h2 = K(E.j_invariant()).global_height_arch()/6
     195    h3 = K(E.b2()/12).global_height_arch()
     196    return n(h1 + h2/2 + h3/2 + c)
     197
     198def Silverman_height_bounds(E):
     199    K = E.base_field()
     200    mu = Faltings_height(E)
     201    lb = -mu-2.14
     202    ub = abs_log_height([K(E.j_invariant()), K(1)])/12 + mu + 1.922
     203    return lb, ub
     204
     205def ReducedQ(E, f, LGen, Elog, C9, C10, Periods, expTors, initQ):
     206    RR=RealField()
     207    w1, w2 = Periods
     208    r = len(LGen)
     209    newQ = initQ
     210    # Repeat LLL reduction until no further reduction is possible
     211    while True:
     212        Q = newQ
     213        S = r*(Q**2)*(expTors**2)
     214        T = 3*r*Q*expTors/sqrt(2.0)
     215        # Create the basis matrix
     216        C = 1
     217        while True:
     218            C *= Q**integer_ceil((r+2)/2)   
     219            precbits = int(RR(C).log(2)+10)
     220            L = copy(zero_matrix(ZZ, r+2))
     221            # Elliptic logarithm should have precision "suitable to" C
     222            # e.g. If C = 10**100, then Re(Elog[i]) should be computed
     223            # correctly to at least 100 decimal places
     224            if precbits > Elog[0].prec():
     225                verbose( "Need higher precision, recompute elliptic logarithm ...")
     226                # Re-compute elliptic logarithm to the right precision
     227                verbose( "precision in bits %i" % precbits)
     228                Elog = [ P.elliptic_logarithm(f, precision = precbits) for P in LGen]
     229                verbose( "Elliptic logarithm recomputed")
     230                w1, w2 = E.period_lattice(f).normalised_basis( prec = precbits)
     231            # Assign all non-zero entries
     232            for i in range(r):
     233                L[i, i] = 1
     234                L[r, i]   = (C*Elog[i].real_part()).trunc()
     235                L[r+1, i] = (C*Elog[i].imag_part()).trunc()
     236            L[r , r ]   = (C*w1.real_part()).trunc()
     237            L[r , r+1 ] = (C*w2.real_part()).trunc()
     238            L[r+1, r]   = (C*w1.imag_part()).trunc()
     239            L[r+1, r+1] = (C*w2.imag_part()).trunc()
     240            # LLL reduction and constants
     241            L = L.transpose()
     242            L = L.LLL()
     243            b1 = L[0] # 1st row of reduced basis
     244            # Norm(b1) = square of Euclidean norm
     245            normb1 = sum([i**2 for i in b1])
     246            lhs = 2**(-r-1)*normb1 - S
     247            if (lhs >= 0) and (sqrt(lhs) > T):
     248                break
     249        newQ = ( RR(C*C9*expTors).log() - RR(sqrt(lhs) - T).log() ) / C10
     250        newQ = integer_floor(sqrt(newQ))
     251        verbose( "After reduction, Q <= %f" % newQ)
     252        if ((Q - newQ) <= 1) or (newQ <= 1):
     253            break
     254    return newQ
     255
     256
     257#// Search for all integral points on elliptic curves over number fields
     258#// within a given bound
     259#// Input:    E = elliptic curve over a number field K
     260#//           L = a sequence of points in the Mordell-Weil basis for E(K)
     261#//           Q = a maximum for the absolute bound on all coefficients
     262#//               in the linear combination of points in L
     263#// Output:  S1 = sequence of all integral points on E(K) modulo [-1]
     264#//          S2 = sequence of tuples representing the points as a
     265#//               linear combination of points in L
     266#// Option: tol = tolerance used for checking integrality of points.
     267#//               (Default = 0 - only exact arithmetic will be performed)
     268
     269# cyclic group iterator 
     270# Returns all elements of the cyclic group, remembering all intermediate steps
     271# id - identity element
     272# gens - generators of the group
     273# mult - orders of the generators
     274# both_signs - whether to return all elements or one per each {element, inverse} pair
     275def cyc_iter(id, gens, mult, both_signs = False):
     276    if len(gens) == 0:
     277        yield id,[]
     278    else:
     279        P = gens[0]
     280        cur = id
     281        if both_signs:
     282            ran = xrange(mult[0])
     283        else:
     284            ran = xrange(mult[0]/2 + 1)
     285        for k in ran:
     286            for rest, coefs in cyc_iter(id, gens[1:], mult[1:], both_signs or k != 0):
     287                yield cur + rest, [k] + coefs
     288            cur += P
     289
     290#generates elements of form a_1G_1 + ... + a_nG_n
     291#where |a_i| <= bound and the first non-zero coefficient is positive
     292def L_points_iter(id, gens, bound, all_zero = True):
     293    if len(gens) == 0:
     294        yield id, []
     295    else:
     296        P = gens[0]
     297        cur = id
     298        for k in xrange(bound+1):
     299            for rest, coefs in L_points_iter(id, gens[1:], bound, all_zero = (all_zero and k == 0)):
     300                yield cur + rest, [k] + coefs 
     301                if k!=0 and not all_zero:
     302                   yield -cur + rest, [-k] + coefs
     303            cur += P
     304
     305def integral_points_with_Q(E, L, Q, tol = 0):
     306    r'''Given an elliptic curve E over a number field, its Mordell-Weil basis L, and a positive integer Q, return the sequence of all integral points modulo [-1] of the form P = q_1*L[1] + ... + q_r*L[r] + T with some torsion point T and |q_i| <= Q, followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points.
     307    '''
     308    assert tol >= 0
     309    CC=ComplexField()
     310    if L == []:
     311        print "Generators must be supplied!"
     312        return None
     313
     314    Tors = E.torsion_subgroup()
     315    expTors = Tors.exponent()
     316    OrdG = Tors.invariants()
     317    Tgens = Tors.gens()
     318    id = E([0,1,0])
     319    total_gens = L + list(Tgens)
     320    verbose( "Generators = %s" % L)
     321    PtsList = []
     322    if tol == 0:
     323        verbose( "Exact arithmetic")
     324        for P, Pcoeff in L_points_iter(id, L, Q):
     325            for T, Tcoeff in cyc_iter(id, Tgens, OrdG, both_signs = P!=id):
     326                R = P + T
     327                if R[0].is_integral() and R[1].is_integral() and R != id:
     328                    Rcoeff = Pcoeff + Tcoeff
     329                    verbose( "%s ---> %s" %(R, Rcoeff))
     330                    PtsList.append(R)
     331        verbose( "*"*45)
     332        return PtsList
     333    # Suggested by John Cremona
     334    # Point search. This is done via arithmetic on complex points on each
     335    # embedding of E. Exact arithmetic will be carried if the resulting
     336    # complex points are "close" to being integral, subject to some tolerance
     337    K = E.base_ring()
     338    r, s = K.signature()
     339    # Set tolerance - This should be larger than 10**(-current precision) to
     340    # avoid missing any integral points. Too large tolerance will slow the
     341    # computation, too small tolerance may lead to missing some integral points.
     342    verbose( "Tolerance = %f" % tol )
     343    if K == QQ:
     344        A = matrix([1])
     345    else:
     346        A = matrix([a.complex_embeddings() for a in K.integral_basis()])
     347        # Note that A is always invertible, so we can take its inverse
     348    B = A.inverse()
     349    point_dict = {}
     350    for emb in K.embeddings(CC):
     351        Ec = EllipticCurve(CC, map(emb,E.a_invariants()))
     352        # need to turn off checking, otherwise the program crashes
     353        Lc = [ Ec.point(map(emb,[P[0],P[1],P[2]]), check = False) for P in L]
     354        Tgensc = [ Ec.point(map(emb,[P[0],P[1],P[2]]), check = False) for P in Tgens]
     355        # Compute P by complex arithmetic for every embedding
     356        id = Ec([0,1,0])
     357        for P, Pcoeff in L_points_iter(id, Lc, Q):
     358            for T, Tcoeff in cyc_iter(id, Tgensc, OrdG, both_signs = P!=id):
     359                R = P + T
     360                if R == id:
     361                    continue
     362                key = tuple(Pcoeff + Tcoeff)
     363                if key in point_dict:
     364                    point_dict[key] += [R]
     365                else:
     366                    point_dict[key] = [R]
     367    integral_pts = []
     368    false = 0
     369    for Pcoeff in point_dict:
     370        # Check if the x-coordinate of P is "close to" being integral
     371        # If so, compute P exactly and check if it is integral skip P otherwise
     372        XofP = vector([Pemb[0] for Pemb in point_dict[Pcoeff]])
     373        # Write x(P) w.r.t. the integral basis of (K)
     374        # Due to the floating arithmetic, some entries in LX may have very tiny
     375        # imaginary part, which can be thought as zero
     376        LX = XofP * B
     377        try:
     378            LX = [abs( i.real_part() - i.real_part().round() ) for i in LX[0]]
     379        except AttributeError:
     380            LX = [abs(i - i.round()) for i in LX[0]]
     381        if len([1 for dx in LX if dx >= tol]) == 0 :
     382        # x-coordinate of P is not integral, skip P
     383            P = sum([Pcoeff[i]*Pt for i, Pt in enumerate(total_gens)])
     384            if P[0].is_integral() and P[1].is_integral():
     385                integral_pts.append(P)
     386                verbose( "%s ---> %s" %(P, Pcoeff))
     387            else:
     388                false += 1
     389    verbose( 'false positives %i'% false)
     390    return integral_pts
     391# Compute all integral points on elliptic curve over a number field
     392# This function simply computes a suitable bound Q, and return
     393# IntegralPoints(E, L, Q    tol = ...).
     394# Input        E = elliptic curve over a number field K
     395#                     L = a sequence of points in the Mordell-Weil basis for E(K)
     396# Output    S1 = sequence of all integral points on E(K) modulo [-1]
     397#                    S2 = sequence of tuples representing the points as a
     398#                             linear combination of points in L
     399# Option tol = tolerance used for checking integrality of points.
     400#                             (Default = 0 - only exact arithmetic will be performed)
     401
     402def calculate_Q(E,L):
     403    from math import pi as Rpi
     404    CC=ComplexField()
     405   
     406    if len(L) == 0:
     407        return 0
     408    K = E.base_ring()
     409    expTors = E.torsion_subgroup().exponent()
     410    r, s = K.signature()
     411    RR = RealField()
     412    pi = RR(Rpi)
     413    b2 = E.b_invariants()[0]
     414    # Global constants (independent to the embedding of E)
     415    gl_con = {}
     416    gl_con['c2'] = C2 = - Silverman_height_bounds(E)[0]
     417    gl_con['c3'] = C3 = c3(L)
     418    gl_con['h_E'] = h = h_E(E)
     419    verbose( "Global constants")
     420    for name, val in gl_con.iteritems():
     421        verbose( '%s = %s'%(name, val))
     422    verbose( "-"*45)
     423    Q = []
     424    # Find the most reduced bound on each embedding of E
     425    # Sage bug, QQ.places() is not implemented
     426    # and K.embeddings gives each complex places twice
     427    if K is QQ:
     428        infplaces = K.embeddings(CC) 
     429    else:
     430        infplaces = K.places()
     431    for i, f in enumerate(infplaces):
     432        # Bug in P.complex_logarithm(QQ.embeddings(CC)[0])
     433        if K is QQ:
     434            emb = None
     435        else:
     436            emb = f
     437        if i < r:
     438            nv = 1
     439            verbose( "Real embedding #%i" % i )
     440        else:
     441            nv = 2
     442            verbose( "Complex embedding #%i" % (i-r))
     443        # Create complex embedding of E
     444        # Embedding of all points in Mordell-Weil basis
     445        # Find complex elliptic logarithm of each embedded point
     446        precbits = integer_floor(100*RR(10).log(2))
     447        Elog = [P.elliptic_logarithm(emb, precision = precbits) for P in L]
     448        Periods = E.period_lattice(emb).normalised_basis();
     449        # Local constants (depending on embedding)
     450        # C9, C10 are used for the upper bound on the linear form in logarithm
     451        loc_con = {}
     452        loc_con['c4'] = C4 = C3 * K.degree() / (nv*(r+s))
     453        loc_con['c5'] = C5 = C2 * K.degree() / (nv*(r+s))
     454        loc_con['c6'] = C6 = compute_c6(E,f)
     455        loc_con['delta'] = delta = 1 + (nv-1)*pi
     456        loc_con['c8'] = C8 = compute_c8(Periods)
     457        loc_con['c7'] = C7 = 8*(delta**2) + (C8**2)*abs(f(b2))/12
     458        loc_con['c9'] = C9 = C7*RR(C5/2).exp()
     459        loc_con['c10'] = C10 = C4/2
     460        loc_con['Q0'] = Q0 = sqrt( ( RR(C6+abs(f(b2))/12).log() + C5) / C4 )
     461
     462        # Constants for David's lower bound on the linear form in logarithm
     463        # Magma and Sage output periods in different order, need to swap w1 and w2
     464        w2, w1 = Periods # N.B. periods are already in "standard" form
     465        loc_con['d7'] = D7 = 3*pi / ((abs(w2)**2) * (w2/w1).imag_part())
     466        loc_con['d8'] = D8 = d8(E, L, Elog, Periods, D7)
     467        loc_con['d9'] = D9 = d9(E, L, Elog, Periods, D7)
     468        loc_con['d10'] = D10 = d10(E, L, Elog, Periods, D7)
     469        for name, val in loc_con.iteritems():
     470            verbose( "{0} = {1}".format(name, val))
     471        # Find the reduced bound for the coefficients in the linear logarithmic form
     472        loginitQ = InitialQ(K.degree(), len(L), Q0, C9, C10, D8, D9, D10, h, expTors)
     473        verbose( "Initial Q <= 10^%f" % loginitQ)
     474        initQ = 10**loginitQ
     475        Q.append(ReducedQ(E, emb, L, Elog, C9, C10, Periods, expTors, initQ))
     476        verbose( "-"*45)
     477    Q = max(Q)
     478    verbose( "Maximum absolute bound on coefficients = %i"% Q)
     479    return Q
     480
     481def integral_points(E, L, tol = 0, both_signs = False):
     482    r'''Given an elliptic curve over a number field and its Mordell-Weil basis, return the sequence of all integral points modulo [-1], followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points. (This function simply computes a suitable Q and call
     483IntegralPoints(E, L, Q tol = ...)
     484'''
     485    assert tol >= 0
     486    if L == []:
     487        print "Generators must be supplied!"
     488        return None
     489
     490    id = E([0,1,0])
     491    Q = calculate_Q(E,L)
     492    int_points = integral_points_with_Q(E, L, Q, tol = tol)
     493    if both_signs:
     494        int_points += [-P for P in int_points]
     495    int_points.sort()
     496    return int_points