Ticket #10973: ell_int_points.py

File ell_int_points.py, 19.3 KB (added by Justin Walker, 12 years ago)

This code goes into "sage/schemes/elliptiic_curves". It works, but needs work.

Line 
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.other import sqrt
9from sage.matrix.constructor import zero_matrix
10from sage.matrix.matrix_space import MatrixSpace
11from sage.misc.all import verbose
12from sage.misc.functional import n
13from sage.misc.misc_c import prod
14from sage.rings.integer_ring import ZZ
15from sage.rings.real_mpfr import RealField
16from sage.rings.complex_field import ComplexField
17from sage.rings.rational_field import QQ
18from sage.rings.arith import integer_ceil, integer_floor
19from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
20
21# This function should be detached and included as part of projective
22#  points over number field
23def abs_log_height(X, gcd_one = True, precision = None):
24    r''' Computes the height of a point in a projective space over field K.
25    It assumes the coordinates have gcd equal to 1.
26    If not use gcd_one flag.
27    '''
28    assert isinstance(X,list)
29    K = X[0].parent()
30    if precision is None:
31        RR = RealField()
32        CC = ComplexField()
33    else:
34        RR = RealField(precision)
35        CC = ComplexField(precision)
36    places = set([])
37    if K == QQ:
38        embs = K.embeddings(RR)
39        Xideal = X
40    else:
41        embs = K.places(precision)
42        # skipping zero as it currently over K breaks Sage
43        Xideal = [K.ideal(x) for x in X if not x.is_zero()]
44    for x in Xideal:
45        places = places.union(x.denominator().prime_factors())
46        if not gcd_one:
47            places = places.union(x.numerator().prime_factors())
48    if K == QQ:
49        non_arch = sum([max([RR(P)**(-x.valuation(P)) for x in X]).log() for P in places])
50    else:
51        non_arch = sum([P.residue_class_degree() *
52                        P.absolute_ramification_index() *
53                        max([x.abs_non_arch(P, precision) for x in X]).log() for P in places])
54    arch = 0
55    r,s = K.signature()
56    for i,f in enumerate(embs):
57        if i<r:
58            arch += max([f(x).abs() for x in X]).log()
59        else:
60            arch += max([f(x).abs()**2 for x in X]).log()
61    return (arch+non_arch)/K.degree()
62   
63def compute_c6(E,emb):
64    CC=ComplexField()
65    CCx=PolynomialRing(CC,'x')
66    x = CCx.gen()
67    #f = x**3-27*emb(E.c4())*x-54*emb(E.c6())
68    f = x**3-emb(E.c4()/48)*x-emb(E.c6()/864)
69    R = f.roots(multiplicities = False)
70    m = max([CC(r).abs() for r in R])
71    return 2*m
72
73def compute_c8(L):
74    # Original code transformed the lattice generators.
75    # Here we assume they are of standard form.
76    CC=ComplexField()
77
78    w1, w2 = L
79    m = max(CC(w1).abs(), CC(w2).abs(), CC(w1+w2).abs())
80    return m
81
82# There is a Sage trak ticket implementing this. In the future the library function can be used and this one removed
83def height_pairing_matrix(points, precision = None):
84    r = len(points)
85    if precision is None:
86        RR = RealField()
87    else:
88        RR = RealField(precision)
89    M = MatrixSpace(RR, r)
90    mat = M()
91    for j in range(r):
92        mat[j,j] = points[j].height(precision = precision)
93    for j in range(r):
94        for k in range(j+1,r):
95            mat[j,k] = ((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2
96            mat[k,j] = mat[j,k]
97    return mat
98
99def c3(L):
100   import warnings
101   with warnings.catch_warnings():
102       warnings.simplefilter("ignore")
103       return min(map(abs,height_pairing_matrix(L).eigenvalues()))
104
105def h_E(E):
106    K = E.base_field()
107    j = E.j_invariant()
108    c4, c6 = E.c_invariants()
109    g2, g3 = c4/12, c6/216
110    return max(abs_log_height([K(1), g2, g3]), abs_log_height([K(1), j]))
111
112def h_m(E, P, ElogEmbedP, D7):
113    RR = RealField()
114    K = E.base_field()
115    return max([P.height(), h_E(E), D7/K.degree()*abs(RR(ElogEmbedP).log())**2])
116   
117def Extra_h_m(E, Periods, D7):
118    D = E.base_field().degree()
119    h = h_E(E)
120    return map(lambda x: max([0, h, D7/D*abs(x)**2]), Periods)
121
122def d8(E, L, Elog, Periods, D7):
123    RR=RealField()
124    C = [RR(1).exp() * h_E(E)]
125    D = E.base_field().degree()
126    for i in range(len(L)):
127        C.append(h_m(E, L[i], Elog[i], D7) / D)
128    C += [t / D for t in Extra_h_m(E, Periods, D7)]
129    return max(C);
130
131def d9(E, L, Elog, Periods, D7):
132    RR=RealField()
133    D = E.base_field().degree()
134    C = []
135    for i in range(len(L)):
136        tmp = RR(1.0).exp() * sqrt(D * h_m(E, L[i], Elog[i], D7) / D7)
137        C.append( tmp / abs(Elog[i]))
138    Ehm = Extra_h_m(E, Periods, D7)
139    C += [RR(1).exp() * sqrt(D*Ehm[i]/D7) / abs(Periods[i]) for i in [0,1]]
140    verbose("d9: C = %s"%(C,))
141    return min(C);
142
143def d10(E, L, Elog, Periods, D7):
144    RR=RealField()
145    D = E.base_field().degree()
146    n = len(L)+2
147    scalar = 2 * 10**(8 + 7*n) * (2/RR(1).exp())**(2*n**2)
148    scalar *= (n+1)**(4*n**2 + 10*n) * D**(2*n + 2)
149    scalar *= (RR(d9(E, L, Elog, Periods, D7)).log())**(-2*n-1)
150    for i in range(len(L)):
151        scalar *= h_m(E, L[i], Elog[i], D7)
152    scalar *= prod(Extra_h_m(E, Periods, D7))   
153    return scalar
154
155def RHS(D, r, C9, C10, D9, D10, h, Q, expTors):
156    RR=RealField()
157    bound = (RR(RR(Q*r*expTors).log()).log() + h + RR(D*D9).log())**(r+3)
158    bound *= D10*(RR(Q*r*expTors).log() + RR(D*D9).log())
159    bound += RR(C9*expTors).log()
160    bound /= C10
161    return bound
162
163def InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors):
164    RR=RealField()
165    minQ = max(Q0, RR(D8).exp())
166    Q = minQ + 1
167    x = integer_ceil(Q.log(10)) # x = log_10(Q)
168    verbose("x: %s, %s"%(type(Q), type(x)))
169    exp_OK = 0   # the exponent that satisfies P.I.
170    exp_fail = 0 # the exponent that fails P.I.
171    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))
172    while 10**(2*x) < RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
173        exp_OK = x # Principal Inequality satisfied
174        x *= 2     # double x, and retry
175    exp_fail = x # x that fails the Principal Inequality
176   
177   
178    # So now x = log_10(Q) must lie between exp_OK and exp_fail
179    # Refine x further using "binary search"
180    while True:
181        x = integer_floor((exp_OK + exp_fail)/2)
182        if 10**(2*x) >= RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors): 
183            exp_fail = x
184        else:
185            exp_OK = x
186        if (exp_fail - exp_OK) <= 1:
187            break
188    return exp_fail # over-estimate
189
190def Faltings_height(E):
191    RR=RealField()
192    K = E.base_field()
193    c = RR(2).log()
194    if E.b2().is_zero():
195        c = 0
196    h1 = abs_log_height([K(E.discriminant()), K(1)])/6
197    h2 = K(E.j_invariant()).global_height_arch()/6
198    h3 = K(E.b2()/12).global_height_arch()
199    return n(h1 + h2/2 + h3/2 + c)
200
201def Silverman_height_bounds(E):
202    K = E.base_field()
203    mu = Faltings_height(E)
204    lb = -mu-2.14
205    ub = abs_log_height([K(E.j_invariant()), K(1)])/12 + mu + 1.922
206    return lb, ub
207
208def ReducedQ(E, f, LGen, Elog, C9, C10, Periods, expTors, initQ):
209    RR=RealField()
210    w1, w2 = Periods
211    r = len(LGen)
212    newQ = initQ
213    # Repeat LLL reduction until no further reduction is possible
214    while True: 
215        Q = newQ
216        S = r*(Q**2)*(expTors**2)
217        T = 3*r*Q*expTors/sqrt(2)
218        # Create the basis matrix
219        C = 1
220        while True: 
221            C *= Q**integer_ceil((r+2)/2)   
222            precbits = int(RR(C).log(2)+10)
223            L = copy(zero_matrix(ZZ, r+2))
224            # Elliptic logarithm should have precision "suitable to" C
225            # e.g. If C = 10**100, then Re(Elog[i]) should be computed
226            # correctly to at least 100 decimal places
227            if precbits > Elog[0].prec(): 
228                verbose( "Need higher precision, recompute elliptic logarithm ...")
229                # Re-compute elliptic logarithm to the right precision
230                verbose( "precision in bits %i" % precbits)
231                Elog = [ P.elliptic_logarithm(f, precision = precbits) for P in LGen]
232                verbose( "Elliptic logarithm recomputed")
233                w1, w2 = E.period_lattice(f).normalised_basis( prec = precbits)
234            # Assign all non-zero entries
235            for i in range(r): 
236                L[i, i] = 1
237                L[r, i]   = (C*Elog[i].real_part()).trunc()
238                L[r+1, i] = (C*Elog[i].imag_part()).trunc()
239            L[r , r ]   = (C*w1.real_part()).trunc()
240            L[r , r+1 ] = (C*w2.real_part()).trunc()
241            L[r+1, r]   = (C*w1.imag_part()).trunc()
242            L[r+1, r+1] = (C*w2.imag_part()).trunc()
243            # LLL reduction and constants
244            L = L.transpose()
245            L = L.LLL()
246            b1 = L[0] # 1st row of reduced basis
247            # Norm(b1) = square of Euclidean norm
248            normb1 = sum([i**2 for i in b1])
249            lhs = 2**(-r-1)*normb1 - S
250            if (lhs >= 0) and (sqrt(lhs) > T):
251                break
252        newQ = ( RR(C*C9*expTors).log() - RR(sqrt(lhs) - T).log() ) / C10
253        newQ = integer_floor(sqrt(newQ))
254        verbose( "After reduction, Q <= %f" % newQ)
255        if ((Q - newQ) <= 1) or (newQ <= 1):
256            break
257    return newQ
258
259
260#// Search for all integral points on elliptic curves over number fields
261#// within a given bound
262#// Input:    E = elliptic curve over a number field K
263#//           L = a sequence of points in the Mordell-Weil basis for E(K)
264#//           Q = a maximum for the absolute bound on all coefficients
265#//               in the linear combination of points in L
266#// Output:  S1 = sequence of all integral points on E(K) modulo [-1]
267#//          S2 = sequence of tuples representing the points as a
268#//               linear combination of points in L
269#// Option: tol = tolerance used for checking integrality of points.
270#//               (Default = 0 - only exact arithmetic will be performed)
271
272# cyclic group iterator 
273# Returns all elements of the cyclic group, remembering all intermediate steps
274# id - identity element
275# gens - generators of the group
276# mult - orders of the generators
277# both_signs - whether to return all elements or one per each {element, inverse} pair
278def cyc_iter(id, gens, mult, both_signs = False):
279    if len(gens) == 0:
280        yield id,[] 
281    else:
282        P = gens[0]
283        cur = id
284        if both_signs:
285            ran = xrange(mult[0])
286        else:
287            ran = xrange(mult[0]/2 + 1)
288        for k in ran: 
289            for rest, coefs in cyc_iter(id, gens[1:], mult[1:], both_signs or k != 0):
290                yield cur + rest, [k] + coefs
291            cur += P
292
293#generates elements of form a_1G_1 + ... + a_nG_n
294#where |a_i| <= bound and the first non-zero coefficient is positive
295def L_points_iter(id, gens, bound, all_zero = True):
296    if len(gens) == 0:
297        yield id, []
298    else:
299        P = gens[0]
300        cur = id
301        for k in xrange(bound+1):
302            for rest, coefs in L_points_iter(id, gens[1:], bound, all_zero = (all_zero and k == 0)):
303                yield cur + rest, [k] + coefs 
304                if k!=0 and not all_zero:
305                   yield -cur + rest, [-k] + coefs
306            cur += P
307
308def integral_points_with_Q(E, L, Q, tol = 0):
309    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.
310    '''
311    assert tol >= 0 
312    CC=ComplexField()
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    id = E([0,1,0]) 
487    Q = calculate_Q(E,L)
488    int_points = integral_points_with_Q(E, L, Q, tol = tol)
489    if both_signs:
490        int_points += [-P for P in int_points]
491    int_points.sort()
492    return int_points