Ticket #3674: 3674-jcremona-integral-points.patch

a/sage/schemes/elliptic_curves/ell_generic.py

@@ -497 +497 @@
-            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 

a/sage/schemes/elliptic_curves/ell_point.py

@@ -54 +54 @@
-#                  http://www.gnu.org/licenses/
-#*****************************************************************************
-
-from math import ceil, floor, sqrt
-from sage.structure.element import AdditiveGroupElement, RingElement
-
-import sage.plot.all as plot
@@ -708 +707 @@
-        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.
@@ -729 +754 @@
-            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
+
+            # For some curves (e.g. 3314b3) the default precision is not enough!
+            while len(real_roots)!=1:
+                precision*=2
+                RR=rings.RealField(precision)
+                real_roots = pol.roots(RR,multiplicities=False) 
+            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(rings.ComplexField(precision),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
+            # For some curves (e.g. 2370i5) the default precision is not enough!
+            while len(real_roots)!=3:
+                precision*=2
+                RR=rings.RealField(precision)
+                real_roots = pol.roots(RR,multiplicities=False) 
+
+            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):

a/sage/schemes/elliptic_curves/ell_rational_field.py

@@ -11 +11 @@
-   -- 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
-"""
-
-#*****************************************************************************
@@ -55 +57 @@
-import sage.databases.cremona
-from   sage.libs.pari.all import pari
-import sage.functions.transcendental as transcendental
-
+
-import sage.libs.mwrank.all as mwrank
-import constructor
-from sage.interfaces.all import gp
@@ -84 +86 @@
-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
@@ -93 +93 @@
-Q = RationalField()         
-C = ComplexField()
-R = RealField()
+Z = IntegerRing()
-IR = rings.RealIntervalField(20)
-
-_MAX_HEIGHT=21
@@ -1505 +1506 @@
-        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() == 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
@@ -2313 +2366 @@
-            self.__torsion_order = self.torsion_subgroup().order()
-            return self.__torsion_order
-
-flag=0
+algorithm="pari"
-        """
-        Returns the torsion subgroup of this elliptic curve.
-
-        INPUT:
-
-
-
+
+
+
+
-
-        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()
@@ -2348 +2405 @@
-        try:
-            return self.__torsion_subgroup
-        except AttributeError:
-flag
+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):
-##         """
@@ -3384 +3498 @@
-        we return v = -1.
-
-        INPUT:
-
+f
-                       first discriminant < -4 that satisfies the Heegner
-                       hypothesis
-            verbose (bool) -- (default: True)
@@ -3564 +3678 @@
-        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.
+            """
+            from sage.misc.all import cartesian_product_iterator
+            from sage.groups.generic import multiples
+            xs=set()
+            N=H_q
+            
+            def use(P):
+                """
+                Helper function to record x-coord of a point if integral.
+                """
+                if not P.is_zero():
+                    xP = P[0]
+                    if xP.is_integral():
+                        xs.add(xP)
+                        
+            def use_t(R):
+                """
+                Helper function to record x-coords of a point +torsion if integral.
+                """
+                for T in tors_points:
+                    use(R+T)
+                        
+            # We use a naive method when the number of possibilities is small:
+
+            if r==1 or N<=10:
+                for P in multiples(mw_base[0],N+1):
+                    use_t(P)
+                return xs
+            
+            # Otherwise is is very very much faster to first compute
+            # the linear combinations over RR, and only compte them as
+            # rational points if they are approximately integral
+
+            def is_approx_integral(P):
+                return (abs(P[0]-P[0].round()))<0.001 and (abs(P[1]-P[1].round()))<0.001
+            
+            RR = RealField() #(100)
+            ER = self.change_ring(RR)
+            Rgens = [ER.lift_x(P[0]) for P in mw_base]
+            for i in range(r):
+                if abs(Rgens[i][1]-mw_base[i][1])>abs((-Rgens[i])[1]-mw_base[i][1]):
+                    Rgens[i] = -Rgens[i]
+            for ni in cartesian_product_iterator([range(-N,N+1) for i in range(r-1)]+[range(N+1)]):
+                RP=sum([ni[i]*Rgens[i] for i in range(r)],ER(0))
+                for T in tors_points:
+                    if is_approx_integral(RP+ER(T)):
+                        P = sum([ni[i]*mw_base[i] for i in range(r)],T)
+                        use(P)
+            return xs
+            
+            # Below this point we keep earlier experimental code which
+            # is slower when either N or r is large
+            
+#           if r==2:
+#               for P in multiples(mw_base[0],N+1):
+#                   for Q in multiples(mw_base[1],N+1):
+#                       for R in set([Q+P,Q-P]):
+#                           use_t(R)
+#               return xs
+#           
+#           if r==3:
+#               for P in multiples(mw_base[0],N+1):
+#                   for Q in multiples(mw_base[1],N+1):
+#                       PQ = [P+Q,P-Q]
+#                       PQ = set(PQ + [-R for R in PQ]) # {\pm P\pm Q}
+#                       for R in multiples(mw_base[2],N+1):
+#                           for S in PQ:
+#                               use_t(R+S)
+#               return xs
+#           
+#           # general rank
+#           mults=[list(multiples(mw_base[i],N+1)) for i in range(r)]
+#           for i in range(r-1):
+#               mults[i] = [-P for P in mults[i] if not P.is_zero()] + mults[i]
+#           for Pi in cartesian_product_iterator(mults):
+#               use_t(sum(Pi,self(0)))
+#           return xs
+#
+#  older, even slower  code:
+#
+#           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: