Ticket #10973: trac_10973

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