# Ticket #3674: sage-trac3674.patch

File sage-trac3674.patch, 27.9 KB (added by cremona, 10 years ago)

Replaces earlier patch -- applies to 3.0.4

• ## sage/schemes/elliptic_curves/ell_point.py

```# HG changeset patch
# User John Cremona <john.cremona@gmail.com>
# Date 1216636935 -3600
# Node ID 60d3d284fd133619eed07a5c5c078549812ef8a1
# Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
trac#3674: patch implementing integral points and elliptic logs

diff -r 91af4c3f6b92 -r 60d3d284fd13 sage/schemes/elliptic_curves/ell_point.py```
 a AUTHORS: #                  http://www.gnu.org/licenses/ #***************************************************************************** from math import ceil, floor, sqrt from sage.structure.element import AdditiveGroupElement, RingElement import sage.plot.all as plot class EllipticCurvePoint_field(AdditiveG h = self.curve().pari_curve().ellheight([self[0], self[1]]) return rings.RR(h) def is_on_identity_component(self): """ Returns True iff this point is on the identity component of its curve (over an ordered field) INPUT: self -- a point on a curve over any ordered field (e.g. QQ) OUTPUT: True iff the point is on the identity component of the identity EXAMPLES: sage: E=EllipticCurve('5077a1') sage: [E.lift_x(x).is_on_identity_component() for x in range(-3,5)] [False, False, False, False, False, True, True, True] """ if self.is_zero():       # trivial case return True E = self.curve() if E.discriminant() < 0: # only one component return True gx =E.division_polynomial(2) gxd = gx.derivative() gxdd = gxd.derivative() return ( gxd(self[0]) > 0 and gxdd(self[0]) > 0) def xy(self): """ Return the x and y coordinates of this point, as a 2-tuple. class EllipticCurvePoint_field(AdditiveG return self[0], self[1] else: return self[0]/self[2], self[1]/self[2] def elliptic_logarithm(self, precision=53): """ Returns the elliptic logarithm of this point on an elliptic curve defined over the reals INPUT: - precision: a positive integer (default 53) setting the number of bits of precision required ALGORITHM: See -[Co2] Cohen H., A Course in Computational Algebraic Number Theory GTM 138, Springer 1996 AUTHORS: - Michael Mardaus (2008-07) } - Tobias Nagel (2008-07)    } original version from [Co2] - John Cremona (2008-07)    revision following eclib code EXAMPLES: sage: E = EllipticCurve('389a') sage: E.discriminant() > 0 True sage: P = E([-1,1]) sage: P.is_on_identity_component () False sage: P.elliptic_logarithm (96) 0.4793482501902193161295330101 + 0.9858688507758241022112038491*I sage: Q=E([3,5]) sage: Q.is_on_identity_component() True sage: Q.elliptic_logarithm (96) 1.931128271542559442488585220 An example with negative discriminant, and a torsion point: sage: E = EllipticCurve('11a1') sage: E.discriminant() < 0 True sage: P = E([16,-61]) sage: P.elliptic_logarithm (96) 0.2538418608559106843377589233 sage: E.period_lattice().basis()[0] / P.elliptic_logarithm (96) 5.000000000000000000000000000 """ RR = rings.RealField(precision) CC = rings.ComplexField(precision) if self.is_zero(): return CC(0) #Initialize E = self.curve() a1 = RR(E.a1()) a2 = RR(E.a2()) a3 = RR(E.a3()) b2 = RR(E.b2()) b4 = RR(E.b4()) disc = E.discriminant() x = RR(self[0]) y = RR(self[1]) pol = E.division_polynomial(2) real_roots = pol.roots(RR,multiplicities=False) if disc < 0: #Connected Case # Here we use formulae equivalent to those in Cohen, but better # behaved when roots are close together try: assert len(real_roots) == 1 except: raise ValueError, ' no or more than one real root despite disc < 0' e1 = real_roots[0] roots = pol.roots(CC,multiplicities=False) roots.remove(e1) e2,e3 = roots zz = (e1-e2).sqrt() # complex beta = (e1-e2).abs() a = 2*zz.abs() b = 2*zz.real(); c = (x-e1+beta)/((x-e1).sqrt()) while (a - b)/a > 0.5**(precision-1): a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2 z = (a/c).arcsin() w = 2*y+a1*x+a3 if w*((x-e1)*(x-e1)-beta*beta) >= 0: z = RR.pi() - z if w>0: z = z + RR.pi() z /= a return z else:                    #Disconnected Case, disc > 0 real_roots.sort() # increasing order real_roots.reverse() # decreasing order e1>e2>e3 try: assert len(real_roots) == 3 except: raise ValueError, ' no or more than one real root despite disc < 0' e1, e2, e3 = real_roots w1, w2 = E.period_lattice().basis(precision) a = (e1 - e3).sqrt() b = (e1 - e2).sqrt() on_egg = (x < e1) # if P is on the "egg", replace it by P+T3 # where T3=(e3,y3) is a 2-torsion point on # the egg coming from w2/2 on the lattice if on_egg: y3 = -(a1*e3+a3)/2 lam = (y-y3)/(x-e3) x3 = lam*(lam+a1)-a2-x-e3 y = lam*(x-x3)-y-a1*x3-a3 x = x3 c = (x - e3).sqrt() while (a-b)/a > 0.5**precision: a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2 z = (a/c).arcsin()/a if (2*y+a1*x+a3) > 0: z = w1 - z if on_egg: z = z + w2/2 return z ##############################  end  ################################ class EllipticCurvePoint_finite_field(EllipticCurvePoint_field): def _magma_init_(self):
• ## sage/schemes/elliptic_curves/ell_rational_field.py

`diff -r 91af4c3f6b92 -r 60d3d284fd13 sage/schemes/elliptic_curves/ell_rational_field.py`
 a AUTHORS: -- Christian Wuthrich (2007): added padic sha computation -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files. -- John Cremona (2008-01) -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points -- John Cremona (2008-07): further work on integral_points """ #***************************************************************************** import sage.databases.cremona import sage.databases.cremona from   sage.libs.pari.all import pari import sage.functions.transcendental as transcendental import math from sage.calculus.calculus import sqrt, floor, ceil import sage.libs.mwrank.all as mwrank import constructor from sage.interfaces.all import gp import ell_tate_curve import ell_tate_curve factor = arith.factor sqrt = math.sqrt exp = math.exp mul = misc.mul next_prime = arith.next_prime kronecker_symbol = arith.kronecker_symbol Q = RationalField() Q = RationalField() C = ComplexField() R = RealField() Z = IntegerRing() IR = rings.RealIntervalField(20) _MAX_HEIGHT=21 class EllipticCurve_rational_field(Ellip self.__regulator[proof] = reg return self.__regulator[proof] def height_pairing_matrix(self, points=None, precision=None): """ Returns the height pairing matrix of the given points on this curve, which must be defined over Q. INPUT: points -- either a list of points, which must be on this curve, or (default) None, in which case self.gens() will be used. precision -- number of bit of precision of result (default: default RealField) TODO: implement variable precision for heights of points EXAMPLES: sage: E = EllipticCurve([0, 0, 1, -1, 0]) sage: E.height_pairing_matrix() [0.0511114082399688] For rank 0 curves, the result is a valid 0x0 matrix: sage: EllipticCurve('11a').height_pairing_matrix() [] sage: E=EllipticCurve('5077a1') sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]]) [  1.36857250535393  -1.30957670708658 -0.634867157837156] [ -1.30957670708658   2.71735939281229   1.09981843056673] [-0.634867157837156   1.09981843056673  0.668205165651928] """ if points is None: points = self.gens() else: for P in points: assert P.curve() is self r = len(points) if precision is None: RR = rings.RealField() else: RR = rings.RealField(precision) M = matrix.MatrixSpace(RR, r) mat = M() for j in range(r): mat[j,j] = points[j].height() for j in range(r): for k in range(j+1,r): mat[j,k]=((points[j]+points[k]).height() - mat[j,j] - mat[k,k])/2 mat[k,j]=mat[j,k] return mat def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False): """ Given a list of rational points on E, compute the saturation class EllipticCurve_rational_field(Ellip OUTPUT: The EllipticCurveTorsionSubgroup instance associated to this elliptic curve. NOTE: To see the torsion points as a list, use torsion_points() EXAMPLES: sage: EllipticCurve('11a').torsion_subgroup() Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C5 associated to the Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field class EllipticCurve_rational_field(Ellip self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag) self.__torsion_order = self.__torsion_subgroup.order() return self.__torsion_subgroup def torsion_points(self, flag=0): """ Returns the torsion points of this elliptic curve as a sorted list. INPUT: flag -- (default: 0)  chooses PARI algorithm: flag = 0: uses Doud algorithm flag = 1: uses Lutz-Nagell algorithm OUTPUT: A list of all the torsion points on this elliptic curve. EXAMPLES: sage: EllipticCurve('11a').torsion_points() [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)] sage: EllipticCurve('37b').torsion_points() [(0 : 1 : 0), (8 : -19 : 1), (8 : 18 : 1)] sage: E=EllipticCurve([-1386747,368636886]) sage: E.torsion_subgroup() Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C8 x C2 associated to the Elliptic Curve defined by y^2  = x^3 - 1386747*x + 368636886 over Rational Field sage: E.torsion_points() [(-1293 : 0 : 1), (-933 : -29160 : 1), (-933 : 29160 : 1), (-285 : -27216 : 1), (-285 : 27216 : 1), (0 : 1 : 0), (147 : -12960 : 1), (147 : 12960 : 1), (282 : 0 : 1), (1011 : 0 : 1), (1227 : -22680 : 1), (1227 : 22680 : 1), (2307 : -97200 : 1), (2307 : 97200 : 1), (8787 : -816480 : 1), (8787 : 816480 : 1)] """ multiples = sage.groups.generic.multiples gens = self.torsion_subgroup(flag).gens() ngens = len(gens) if ngens == 0: return [self(0)] pts = list(multiples(gens[0],gens[0].order())) if ngens == 1: pts.sort() return pts pts = [P+Q for P in pts for Q in multiples(gens[1],gens[1].order())] pts.sort() return pts ## def newform_eval(self, z, prec): ##         """ class EllipticCurve_rational_field(Ellip Eq = ell_tate_curve.TateCurve(self,p) self._tate_curve[p] = Eq return Eq def height(self, precision=53): """Returns real height of this elliptic curve This is used in integral_points() INPUT: precision (integer) bit precision of the field of real numbers in which the result will lie (default: 53 as in RealField()) EXAMPLES: sage: E=EllipticCurve('5077a1') sage: E.height() 17.4513334798896 sage: E.height(100) 17.451333479889612702508579399 sage: E=EllipticCurve([0,0,0,0,1]) sage: E.height() 1.38629436111989 sage: E=EllipticCurve([0,0,0,1,0]) sage: E.height() 7.45471994936400 """ #Computes height of curve 'self' according to Co1 R = RealField(precision) c4 = self.c4() c6 = self.c6() j = self.j_invariant() log_g2 = R((c4/12)).abs().log() log_g3 = R((c6/216)).abs().log() if j == 0: h_j = R(1) else: h_j = max(log(R(abs(j.numerator()))), log(R(abs(j.denominator())))) if (self.c4() != 0) and (self.c6() != 0): h_gs = max(1, log_g2, log_g3) elif c4 == 0: if c6 == 0: return max(1,h_j) h_gs = max(1, log_g3) else: h_gs = max(1, log_g2) return max(R(1),h_j, h_gs) def integral_points(self, mw_base='auto', both_signs=False): """ Computes all integral points (up to sign) on the elliptic curve E which has Mordell-Weil basis mw_base. INPUT: self -- EllipticCurve_Rational_Field mw_base -- list of EllipticCurvePoint generating the Mordell-Weil group of E (default: 'auto' - calls self.gens()) both_signs -- True/False (default False): if True the output contains both P and -P, otherwise only one of ecah pair. OUTPUT: set of all integral points on E (up to sign unless both_signs is True) HINTS: - The complexity increases exponentially in the rank of curve E. - It can help if you try another Mordell-Weil base, because the computation time depends on this, too. EXAMPLES: A curve of rank 3 with no torsion points sage: E=EllipticCurve([0,0,1,-7,6]) sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6)) sage: a=E.integral_points([P1,P2,P3]); a [(-3 : 0 : 1), (-2 : 3 : 1), (-1 : 3 : 1), (0 : 2 : 1), (1 : 0 : 1), (2 : 0 : 1), (3 : 3 : 1), (4 : 6 : 1), (8 : 21 : 1), (11 : 35 : 1), (14 : 51 : 1), (21 : 95 : 1), (37 : 224 : 1), (52 : 374 : 1), (93 : 896 : 1), (342 : 6324 : 1), (406 : 8180 : 1), (816 : 23309 : 1)] It is also possible to not specify mw_base and tors_base but then the Mordell-Weil basis must be found by other functions which will takes longer sage: E=EllipticCurve([0,0,1,-7,6]) sage: a=E.integral_points(both_signs=True); a [(-3 : -1 : 1), (-3 : 0 : 1), (-2 : -4 : 1), (-2 : 3 : 1), (-1 : -4 : 1), (-1 : 3 : 1), (0 : -3 : 1), (0 : 2 : 1), (1 : -1 : 1), (1 : 0 : 1), (2 : -1 : 1), (2 : 0 : 1), (3 : -4 : 1), (3 : 3 : 1), (4 : -7 : 1), (4 : 6 : 1), (8 : -22 : 1), (8 : 21 : 1), (11 : -36 : 1), (11 : 35 : 1), (14 : -52 : 1), (14 : 51 : 1), (21 : -96 : 1), (21 : 95 : 1), (37 : -225 : 1), (37 : 224 : 1), (52 : -375 : 1), (52 : 374 : 1), (93 : -897 : 1), (93 : 896 : 1), (342 : -6325 : 1), (342 : 6324 : 1), (406 : -8181 : 1), (406 : 8180 : 1), (816 : -23310 : 1), (816 : 23309 : 1)] An example with negative discriminant: sage: EllipticCurve('900d1').integral_points() [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)] Another example with rank 5 and no torsion points sage: E=EllipticCurve([-879984,319138704]) sage: P1=E.point((540,1188)); P2=E.point((576,1836)) sage: P3=E.point((468,3132)); P4=E.point((612,3132)) sage: P5=E.point((432,4428)) sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (400s!) 54 NOTES: - This function uses the algorithm given in [Co1] REFERENCES: -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations GTM 239, Springer 2007 AUTHORS: - Michael Mardaus (2008-07) - Tobias Nagel (2008-07) - John Cremona (2008-07) """ ##################################################################### # INPUT CHECK ####################################################### try: assert self.is_integral() except: raise ValueError, "integral_points() can only be called on an integral model" if mw_base=='auto': mw_base = self.gens() r = len(mw_base) else: try: r = len(mw_base) except TypeError: raise TypeError, 'mw_base must be a list' for P in mw_base: try: assert P.curve() is self except: raise ValueError, "points are not on the correct curve" tors_base = self.torsion_subgroup().gens() len_tors = len(tors_base) tors_points = self.torsion_points() # END INPUT-CHECK#################################################### ##################################################################### ##################################################################### # INTERNAL FUNCTIONS ################################################ ############################## begin ################################ def point_preprocessing(list): #Transforms mw_base so that at most one point is on the #compact component of the curve Q = [] mem = -1 for i in range(0,len(list)): if not list[i].is_on_identity_component(): # i.e. is on "egg" if mem == -1: mem = i else: Q.append(list[i]+list[mem]) mem = i else: Q.append(list[i]) if mem != -1: #add last point, which is not in egg, to Q Q.append(2*list[mem]) return Q ##############################  end  ################################ ############################## begin ################################ def modified_height(i):#computes modified height if base point i return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2) ##############################  end  ################################ ############################## begin ################################ def search_points(xmin,xmax): """Returns a list of integral points on self with x in the interval """ if both_signs: return sum([self.lift_x(x,all=True) for x in range(xmin,xmax)],[]) else: return sum([self.lift_x(x,all=True)[1:] for x in range(xmin,xmax)],[]) ##############################  end  ################################ ############################## begin ################################ def search_remaining_points(): """Returns list of integral points on curve E written as linear combination of n times the mordell-weil base points and torsion points (n is bounded by H_q, which will be computed at runtime) """ OP=self(0) pts=[] for i in range(1,((2*H_q+1)**r)/2): koeffs = Z(i).digits(base=2*H_q+1) koeffs = [0]*(r-len(koeffs)) + koeffs P = self(0) for index in range(r): P = P + (koeffs[index]-H_q)*mw_base[index] for Q in tors_points: # P + torsion points (includes 0) P = P + Q if P!=OP and P[0].is_integral(): pts.append(P) if both_signs: Q = -P if not Q in pts: pts.append(Q) return pts ##############################  end  ################################# # END Internal functions ############################################# ###################################################################### if (r == 0): pts = self.torsion_points() pts = [P for P in pts if P[0].is_integral()] pts.sort() return pts e = R(1).exp() pi = R(constants.pi) g2 = self.c4()/12 g3 = self.c6()/216 disc = self.discriminant() j = self.j_invariant() b2 = self.b2() Qx = rings.PolynomialRing(RationalField(),'x') pol = Qx([-self.c6()/216,-self.c4()/12,0,4]) if disc > 0: # two real component -> 3 roots in RR e1,e2,e3 = pol.roots(R,multiplicities=False) mw_base = point_preprocessing(mw_base) #at most one point in E^{egg} elif disc < 0: # one real component => 1 root in RR (=: e3), # 2 roots in C (e1,e2) roots = pol.roots(C,multiplicities=False) e3 = pol.roots(R,multiplicities=False)[0] roots.remove(e3) e1,e2 = roots # Algorithm presented in [Co1] h_E = self.height() w1, w2 = self.period_lattice().basis() mu = R(disc).abs().log() if j!=0: mu += max(R(1),R(j).abs().log()) if b2!=0: mu += max(R(1),R(b2).abs().log()) mu += log(R(2)) c1 = (mu + 2.14).exp() M = self.height_pairing_matrix(mw_base) c2 = min(M.charpoly ().roots(multiplicities=False)) c3 = (w1**2)*R(b2).abs()/48 + 8 c5 = (c1*c3).sqrt() c7 = abs((3*pi)/((w1**2) * (w1/w2).imag())) mw_base_log = [] #contains \Phi(Q_i) mod_h_list = []  #contains h_m(Q_i) c9_help_list = [] for i in range(0,r): mw_base_log.append(mw_base[i].elliptic_logarithm().abs()) mod_h_list.append(modified_height(i)) c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i]) c8 = max(e*h_E,max(mod_h_list)) c9 = e/c7.sqrt() * min(c9_help_list) n=r+1 c10 = R(2 * 10**(8+7*n) * R((2/e)**(2 * n**2)) * (n+1)**(4 * n**2 + 10 * n) * log(c9)**(-2*n - 1) * misc.prod(mod_h_list)) top = 128 #arbitrary first upper bound bottom = 0 log_c9=log(c9); log_c5=log(c5) log_r_top = log(R(r*(10**top))) while R(c10*(log_r_top+log_c9)*(log(log_r_top)+h_E+log_c9)**(n+1)) > R(c2/2 * (10**top)**2 - log_c5): #initial bound 'top' too small, upshift of search interval bottom = top top = 2*top while top >= bottom: #binary-search like search for fitting exponent (bound) bound = floor(bottom + (top - bottom)/2) log_r_bound = log(R(r*(10**bound))) if R(c10*(log_r_bound+log_c9)*(log(log_r_bound)+h_E+log_c9)**(n+1)) > R(c2/2 * (10**bound)**2 - log_c5): bottom = bound + 1 else: top = bound - 1 H_q = R(10)**bound break_cond = 0 #at least one reduction step #reduction via LLL while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10% c = R((H_q**n)*10)  #c has to be greater than H_q^n M = matrix.MatrixSpace(Z,n) m = M.identity_matrix() for i in range(r): m[i, r] = R(c*mw_base_log[i]).round() m[r,r] = R(c*w1).round() #LLL - implemented in sage - operates on rows not on columns m_LLL = m.LLL() m_gram = m_LLL.gram_schmidt()[0] b1_norm = R(m_LLL.row(0).norm()) #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0 c1_LLL = -1 for i in range(n): tmp = R(b1_norm/(m_gram.row(i).norm())) if tmp > c1_LLL: c1_LLL = tmp if c1_LLL < 0: raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base' d_L_0 = R(b1_norm**2 / c1_LLL) #Reducing of upper bound Q = r * H_q**2 T = (1 + (3/2*r*H_q))/2 if d_L_0 < R(T**2+Q): d_L_0 = 10*(T**2*Q) low_bound = R((sqrt(d_L_0 - Q) - T)/c) #new bound according to low_bound and upper bound [c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3 H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2))) H_q_new = ceil(H_q_new) break_cond = R(H_q_new/H_q) H_q = H_q_new #END LLL-Reduction loop int_points = [] ##for a more detailed view how the function works uncomment ##all print - statements below if disc > 0: ##Points in egg have e2>=x>=e1 int_points += search_points(ceil(e1), floor(e2)+1) #print 'Points on compact component \n',list(int_points) ##Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3 int_points += search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3)))) #print 'Points on compact component and non-compact component part 1 \n', list(int_points) #print 'starting search of remainig points using bound ',H_q int_points += search_remaining_points() #print 'Integral points:\n',list(int_points) #print 'Total amount of points:',len(int_points) if both_signs: int_points = list(set(int_points)) # remove duplicates else: xlist = set([P[0] for P in int_points]) int_points = [self.lift_x(x) for x in xlist] int_points.sort() return int_points def cremona_curves(conductors):