# Ticket #3674: sage-trac3674a.patch

File sage-trac3674a.patch, 28.3 KB (added by cremona, 10 years ago)

Based on 3.0.4

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

```# HG changeset patch
# User John Cremona <john.cremona@gmail.com>
# Date 1216537839 -3600
# Node ID ada4f959e6e5df25190883d9220ed7828be7f2f8
# Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
#3674: preliminary version

diff -r 91af4c3f6b92 -r ada4f959e6e5 sage/schemes/elliptic_curves/ell_point.py```
 a 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(): return True E = self.curve() if E.discriminant() < 0: return False 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 """ 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]) real_roots = E.division_polynomial(2).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] 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: 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 = pi - z if w>0: z = z + 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 ada4f959e6e5 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 """ #***************************************************************************** from sage.libs.pari.all import pari from   sage.libs.pari.all import pari import sage.functions.transcendental as transcendental import math from sage.calculus.calculus import sqrt, arcsin, 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 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, field=RealField()): """Returns real height of this elliptic curve This is used in integral_points() INPUT: field -- a field of real numbers in which the result will lie (default: RealField) EXAMPLES: sage: E=EllipticCurve('5077a1') sage: E.height() 17.4513334798896 sage: E.height(RealField(100)) 17.451333479889612702508579399 sage: E.height(RDF) 17.4513334799 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 = field 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'): """ Computes all integral points 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()) OUTPUT: set of all integral points on E 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 : -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)] 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(); 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)] 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!) 108 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 log_plus(x): x = R(x) try: return max(R(1), log(R(x).abs())) except: ## if log(|x|) fails, i.e. x is 0 return R(1) ##############################  end  ################################ ############################## begin ################################ def extract_realroots(list, eps = 1e-8): #function is used to extract real roots from result of #function roots() listed in list #at the same time cutting off impreciseness of roots() #RETURN: sorted list of roots in RR tmp = [] for z,e in list: im = z.imag() re = z.real() if ((im < eps) and (im > -eps)): #imaginary part in eps if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included) tmp.append(0) else: tmp.append(re) tmp.sort() return tmp ##############################  end  ################################ ############################## begin ################################ def extract_roots(list, eps = 1e-8): #function is used to extract roots from result of #function roots() listed in list #cutting off impreciseness of roots() #RETURN: sorted list of roots in C tmp = [] for i in range(0, len(list)): im = list[i][0].imag() re = list[i][0].real() if ((im < eps) and (im > -eps)): if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included) tmp.append(0 + 0*I) else: tmp.append(re + 0*I) elif ((re < eps) and (re > -eps)): tmp.append(im*I) else: #no impreciseness tmp.append(list[i][0]) return tmp ##############################  end  ################################ ############################## 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 """ return sum([self.lift_x(x,all=True) 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) Q = -P if Q!=P: 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 = exp(1) I = constants.I pi = R(constants.pi) g2 = self.c4()/12 g3 = self.c6()/216 disc = self.discriminant() rdisc = R(disc) j = self.j_invariant() b2 = self.b2() rb2 = R(b2) Qx = rings.PolynomialRing(RationalField(),'x') pol = Qx([-self.c6()/216,-self.c4()/12,0,4]) rootsR = pol.roots(R,multiplicities=False) rootsC = pol.roots(C,multiplicities=False) if disc > 0: # two real component -> 3 roots in RR e1,e2,e3 = rootsR 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 = extract_roots(tmp) if roots[0].imag() == 0: e3 = roots[0].real() e2 = roots[1] e1 = roots[2] elif roots[1].imag() == 0: e3 = roots[1].real() e2 = roots[0] e1 = roots[2] else: e3 = roots[2].real() e2 = roots[1] e1 = roots[0] # Algorithm presented in [Co1] h_E = self.height() w1, w2 = self.period_lattice().basis() mu = (log(rdisc.abs())  + log_plus(j))/6 + log_plus(rb2/12) if b2!=0: mu=mu+log(2.0) c1 = R(exp(mu + 2.14)) M = self.height_pairing_matrix(mw_base) c2 = min(M.charpoly ().roots(multiplicities=False)) c3 = (w1**2)*abs(rb2)/48 + 8 c5 = R(sqrt(c1*c3)) 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/sqrt(c7) * 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) int_points = list(set(int_points)) # remove duplicates int_points.sort() return int_points def cremona_curves(conductors):