Ticket #10973: intpts.sage

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