Ticket #10973: trac_10973.2

File trac_10973.2, 20.9 KB (added by justin, 11 years ago)

Updated patch with code to reject requests w/o generators; replaces previous patch of same name.

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