Ticket #3674: sage-trac3674.patch

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
+            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):

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
@@ -1490 +1491 @@
-        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
@@ -2310 +2363 @@
-        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
@@ -2336 +2392 @@
-            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):
-##         """
@@ -3549 +3657 @@
-        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):