# Ticket #3674: sage-trac3674new.patch

File sage-trac3674new.patch, 35.7 KB (added by cremona, 10 years ago)

Replaces ALL above patches

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

```# HG changeset patch
# User John Cremona <john.cremona@gmail.com>
# Date 1218142976 -3600
# Node ID 00987cd79c91861457154b22b1516357b55cb158
# Parent  5e2426d277ebfd4f7010fb71fee7d5ab3f3574b2
#3674: patch implementing integral points over Q and related functions

diff -r 5e2426d277eb -r 00987cd79c91 sage/schemes/elliptic_curves/ell_generic.py```
 a class EllipticCurve_generic(plane_curve. args = tuple(args[0]) return plane_curve.ProjectiveCurve_generic.__call__(self, *args, **kwds) def is_x_coord(self, x): """ Returns whether x is the x-coordinate of a point on this curve. See also lift_x() to find the point(s) with a given x_coordinate.  This function may be useful in cases where testing an element of tha base field for being a square is faster than finding its square root. EXAMPLES: sage: E = EllipticCurve('37a'); E Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field sage: E.is_x_coord(1) True sage: E.is_x_coord(2) True There are no rational points with x-cordinate 3. sage: E.is_x_coord(3) False However, there are such points in \$E(\R)\$: sage: E.change_ring(RR).is_x_coord(3) True And of course it always works in \$E(\C)\$: sage: E.change_ring(RR).is_x_coord(-3) False sage: E.change_ring(CC).is_x_coord(-3) True AUTHOR: John Cremona, 2008-08-07 (adapted from lift_x()) TEST: sage: E=EllipticCurve('5077a1') sage: [x for x in srange(-10,10) if E.is_x_coord (x)] [-3, -2, -1, 0, 1, 2, 3, 4, 8] sage: F=GF(32,'a') sage: E=EllipticCurve(F,[1,0,0,0,1]) sage: set([P[0] for P in E.points() if P!=E(0)]) == set([x for x in F if E.is_x_coord(x)]) True """ K = self.base_ring() try: x = K(x) except TypeError: raise TypeError, 'x must be coercible into the base ring of the curve' a1, a2, a3, a4, a6 = self.ainvs() fx = ((x + a2) * x + a4) * x + a6 if a1.is_zero() and a3.is_zero(): return fx.is_square() b = (a1*x + a3) if K.characteristic() == 2: R = PolynomialRing(K, 'y') F = R([-fx,b,1]) return len(F.roots())>0 D = b*b + 4*fx return D.is_square() def lift_x(self, x, all=False): """ Given the x-coordinate of a point on the curve, use the defining
• ## sage/schemes/elliptic_curves/ell_point.py

`diff -r 5e2426d277eb -r 00987cd79c91 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-1): 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 5e2426d277eb -r 00987cd79c91 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 self.__torsion_order = self.torsion_subgroup().order() return self.__torsion_order def torsion_subgroup(self, flag=0): def torsion_subgroup(self, algorithm="pari"): """ Returns the torsion subgroup of this elliptic curve. INPUT: flag -- (default: 0)  chooses PARI algorithm: flag = 0: uses Doud algorithm flag = 1: uses Lutz-Nagell algorithm algorithm -- string: "pari"  -- (default) use the pari library "doud" -- use Doud's algorithm "lutz_nagell" -- use the Lutz-Nagell theorem 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() class EllipticCurve_rational_field(Ellip try: return self.__torsion_subgroup except AttributeError: self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag) self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, algorithm) self.__torsion_order = self.__torsion_subgroup.order() return self.__torsion_subgroup def torsion_points(self, algorithm="pari"): """ Returns the torsion points of this elliptic curve as a sorted list. INPUT: algorithm -- string: "pari"  -- (default) use the pari library "doud" -- use Doud's algorithm "lutz_nagell" -- use the Lutz-Nagell theorem 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: T=E.torsion_subgroup(); T 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: T == E.torsion_subgroup(algorithm="doud") True sage: T == E.torsion_subgroup(algorithm="lutz_nagell") True 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(algorithm).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 we return v = -1. INPUT: D (int) -- (deault: 0) Heegner discriminant; if 0, use the D (int) -- (default: 0) Heegner discriminant; if 0, use the first discriminant < -4 that satisfies the Heegner hypothesis verbose (bool) -- (default: True) 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 the real height of this elliptic curve. This is used in integral_points() INPUT: precision -- (integer: default 53) bit precision of the field of real numbers in which the result will lie 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 """ 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 antilogarithm(self, z, precision=100): """ Returns the rational point (if any) associated to this complex number; the inverse of the elliptic logarithm function. INPUT: z -- a complex number representing an element of CC/L where L is the period lattice of the elliptic curve precision -- an integer specifying the precision (in bits) which will be used for the computation OUTPUT: The rational point which is the image of z under the Weierstrass parametrization, if it exists and can be determined from z with default precision. NOTE: This uses the function ellztopoint from the pari library TODO: Extend the wrapping of ellztopoint() to allow passing of the precision parameter. """ E_pari = self.pari_curve(prec) try: coords = E_pari.ellztopoint(z) if len(coords) == 1: return self(0) return self([CC(xy).real().simplest_rational() for xy in coords]) except TypeError: raise ValueError, "No rational point computable from z" def integral_points(self, mw_base='auto', both_signs=False, verbose=False): """ Computes all integral points (up to sign) on this elliptic curve. INPUT: 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 each pair. verbose -- True/False (default False): if True, some details of the computation are output OUTPUT: A sorted list of all the integral points on E (up to sign unless both_signs is True) NOTES: The complexity increases exponentially in the rank of curve E.  The computation time (but not the output!) depends on the Mordell-Weil basis.  If mw_base is given but it not a basis for the Mordell-Weil group (modulo torsion), integral points which are not in the subgroup generated by the given points will almost certainly not be listed. 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)] sage: a = E.integral_points([P1,P2,P3], verbose=True) Using mw_basis  [(2 : 0 : 1), (4 : 6 : 1), (114/49 : -720/343 : 1)] e1,e2,e3:  -3.01243037259331 1.06582054769620 1.94660982489710 Minimal eigenvalue of height pairing matrix:  0.472730555831538 x-coords of points on compact component with  -3 <=x<= 1 set([0, -1, -3, -2, 1]) x-coords of points on non-compact component with  2 <=x<= 6 set([2, 3, 4]) starting search of remaining points using coefficient bound  6 x-coords of extra integral points: set([2, 3, 4, 37, 406, 8, 11, 14, 816, 52, 21, 342, 93]) Total number of integral points: 18 It is also possible to not specify mw_base, but then the Mordell-Weil basis must be computed;  this may take much 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 ####################################################### if not self.is_integral(): 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' if not all([P.curve() is self for P in mw_base]): raise ValueError, "points are not on the correct curve" 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 integral_x_coords_in_interval(xmin,xmax): """ Returns the set of integers x with xmin<=x<=xmax which are x-coordinates of points on this curve. """ return set([x for x  in range(xmin,xmax) if self.is_x_coord(x)]) ##############################  end  ################################ ############################## begin ################################ def integral_points_with_bounded_mw_ceoffs(): r""" Returns the set of integers x which are x-coordinates of points on the curve which are linear combinations of the generators (basis and torsion points) with coefficients bounded by \$H_q\$.  The bound \$H_q\$ will be computed at runtime. """ xs=set() for i in range(ceil(((2*H_q+1)**r)/2)): koeffs = Z(i).digits(base=2*H_q+1) koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0)) for Q in tors_points: # P + torsion points (includes 0) tmp = P + Q if not tmp.is_zero(): x = tmp[0] if x.is_integral(): xs.add(x) return xs ##############################  end  ################################# # END Internal functions ############################################# ###################################################################### if (r == 0): int_points = [P for P in tors_points if not P.is_zero()] int_points = [P for P in int_points if P[0].is_integral()] if not both_signs: xlist = set([P[0] for P in int_points]) int_points = [self.lift_x(x) for x in xlist] int_points.sort() if verbose: print 'Total number of integral points:',len(int_points) return int_points 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 #on curve 897e4, only one root is found with default precision! RR = R prec = RR.precision() ei = pol.roots(RR,multiplicities=False) while len(ei)<3: prec*=2 RR=RealField(prec) ei = pol.roots(RR,multiplicities=False) e1,e2,e3 = ei if r >= 2: #preprocessing of mw_base only necessary if rank > 1 mw_base = point_preprocessing(mw_base) #at most one point in #E^{egg}, saved in P_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 e = R(1).exp() pi = R(constants.pi) if verbose: print "Using mw_basis ",mw_base print "e1,e2,e3: ",e1,e2,e3 # Algorithm presented in [Co1] h_E = self.height() w1, w2 = self.period_lattice().basis() mu = R(disc).abs().log() / 6 if j!=0: mu += max(R(1),R(j).abs().log()) / 6 if b2!=0: mu += max(R(1),R(b2).abs().log()) mu += log(R(2)) else: mu += 1 c1 = (mu + 2.14).exp() M = self.height_pairing_matrix(mw_base) c2 = min(M.charpoly ().roots(multiplicities=False)) if verbose: print "Minimal eigenvalue of height pairing matrix: ", c2 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))) #        if verbose: #            print "[bottom,top] = ",[bottom,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) #            if verbose: #                print "[bottom,top] = ",[bottom,top] 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 if low_bound != 0: H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2))) H_q_new = ceil(H_q_new) if H_q_new == 1: break_cond = 1 # stops reduction else: break_cond = R(H_q_new/H_q) H_q = H_q_new else: break_cond = 1 # stops reduction, so using last H_q > 0 #END LLL-Reduction loop b2_12 = b2/12 if disc > 0: ##Points in egg have X(P) between e1 and e2 [X(P)=x(P)+b2/12]: x_int_points = integral_x_coords_in_interval(ceil(e1-b2_12), floor(e2-b2_12)+1) if verbose: print 'x-coords of points on compact component with ',ceil(e1-b2_12),'<=x<=',floor(e2-b2_12) print x_int_points else: x_int_points = set() ##Points in noncompact component with X(P)< 2*max(|e1|,|e2|,|e3|) , espec. X(P)>=e3 x0 = ceil(e3-b2_12) x1 = ceil(2*max(abs(e1),abs(e2),abs(e3)) - b2_12) x_int_points2 = integral_x_coords_in_interval(x0, x1) x_int_points = x_int_points.union(x_int_points2) if verbose: print 'x-coords of points on non-compact component with ',x0,'<=x<=',x1-1 print x_int_points2 if verbose: print 'starting search of remaining points using coefficient bound ',H_q x_int_points3 = integral_points_with_bounded_mw_ceoffs() x_int_points = x_int_points.union(x_int_points3) if verbose: print 'x-coords of extra integral points:' print x_int_points3 # Now we have all the x-coordinates of integral points, and we # construct the points, depending on the parameter both_signs: if both_signs: int_points = sum([self.lift_x(x,all=True) for x in x_int_points],[]) else: int_points = [self.lift_x(x) for x in x_int_points] int_points.sort() if verbose: print 'Total number of integral points:',len(int_points) return int_points def cremona_curves(conductors):