Ticket #3674: sage-trac3674.patch

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

Replaces earlier patch -- applies to 3.0.4

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1216636935 -3600
    # Node ID 60d3d284fd133619eed07a5c5c078549812ef8a1
    # Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
    trac#3674: patch implementing integral points and elliptic logs
    
    diff -r 91af4c3f6b92 -r 60d3d284fd13 sage/schemes/elliptic_curves/ell_point.py
    a b AUTHORS: 
    5454#                  http://www.gnu.org/licenses/
    5555#*****************************************************************************
    5656
    57 from math import ceil, floor, sqrt
    5857from sage.structure.element import AdditiveGroupElement, RingElement
    5958
    6059import sage.plot.all as plot
    class EllipticCurvePoint_field(AdditiveG 
    708707        h = self.curve().pari_curve().ellheight([self[0], self[1]])
    709708        return rings.RR(h)
    710709   
     710    def is_on_identity_component(self):
     711        """
     712        Returns True iff this point is on the identity component of
     713        its curve (over an ordered field)
     714
     715        INPUT:
     716            self -- a point on a curve over any ordered field (e.g. QQ)
     717
     718        OUTPUT:
     719            True iff the point is on the identity component of the identity
     720
     721        EXAMPLES:
     722            sage: E=EllipticCurve('5077a1')     
     723            sage: [E.lift_x(x).is_on_identity_component() for x in range(-3,5)]
     724            [False, False, False, False, False, True, True, True]
     725        """
     726        if self.is_zero():       # trivial case
     727            return True
     728        E = self.curve()
     729        if E.discriminant() < 0: # only one component
     730            return True
     731        gx =E.division_polynomial(2)
     732        gxd = gx.derivative()
     733        gxdd = gxd.derivative()
     734        return ( gxd(self[0]) > 0 and gxdd(self[0]) > 0)
     735
    711736    def xy(self):
    712737        """
    713738        Return the x and y coordinates of this point, as a 2-tuple.
    class EllipticCurvePoint_field(AdditiveG 
    729754            return self[0], self[1]
    730755        else:
    731756            return self[0]/self[2], self[1]/self[2]
     757
     758    def elliptic_logarithm(self, precision=53):
     759        """ Returns the elliptic logarithm of this point on an
     760            elliptic curve defined over the reals
     761
     762        INPUT: - precision: a positive integer (default 53) setting
     763        the number of bits of precision required
     764
     765        ALGORITHM: See -[Co2] Cohen H., A Course in Computational
     766                        Algebraic Number Theory GTM 138, Springer 1996
     767
     768        AUTHORS:
     769            - Michael Mardaus (2008-07) }
     770            - Tobias Nagel (2008-07)    } original version from [Co2]
     771            - John Cremona (2008-07)    revision following eclib code
     772       
     773        EXAMPLES:
     774            sage: E = EllipticCurve('389a')
     775            sage: E.discriminant() > 0
     776            True
     777            sage: P = E([-1,1])
     778            sage: P.is_on_identity_component ()
     779            False
     780            sage: P.elliptic_logarithm (96)
     781            0.4793482501902193161295330101 + 0.9858688507758241022112038491*I
     782            sage: Q=E([3,5])
     783            sage: Q.is_on_identity_component()     
     784            True
     785            sage: Q.elliptic_logarithm (96)                 
     786            1.931128271542559442488585220
     787
     788            An example with negative discriminant, and a torsion point:
     789            sage: E = EllipticCurve('11a1')
     790            sage: E.discriminant() < 0
     791            True
     792            sage: P = E([16,-61])
     793            sage: P.elliptic_logarithm (96)                               
     794            0.2538418608559106843377589233
     795            sage: E.period_lattice().basis()[0] / P.elliptic_logarithm (96)
     796            5.000000000000000000000000000
     797           
     798        """
     799        RR = rings.RealField(precision)
     800        CC = rings.ComplexField(precision)
     801        if self.is_zero():
     802            return CC(0)
     803
     804        #Initialize
     805
     806        E = self.curve()
     807        a1 = RR(E.a1())
     808        a2 = RR(E.a2())
     809        a3 = RR(E.a3())
     810        b2 = RR(E.b2())
     811        b4 = RR(E.b4())
     812        disc = E.discriminant()
     813        x = RR(self[0])
     814        y = RR(self[1])
     815        pol = E.division_polynomial(2)
     816        real_roots = pol.roots(RR,multiplicities=False)
     817
     818        if disc < 0: #Connected Case
     819            # Here we use formulae equivalent to those in Cohen, but better
     820            # behaved when roots are close together
     821            try:
     822                assert len(real_roots) == 1
     823            except:
     824                raise ValueError, ' no or more than one real root despite disc < 0'
     825            e1 = real_roots[0]
     826            roots = pol.roots(CC,multiplicities=False)
     827            roots.remove(e1)
     828            e2,e3 = roots
     829            zz = (e1-e2).sqrt() # complex
     830            beta = (e1-e2).abs()
     831            a = 2*zz.abs()
     832            b = 2*zz.real();
     833            c = (x-e1+beta)/((x-e1).sqrt())
     834            while (a - b)/a > 0.5**(precision-1):
     835                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     836            z = (a/c).arcsin()
     837            w = 2*y+a1*x+a3           
     838            if w*((x-e1)*(x-e1)-beta*beta) >= 0:
     839                z = RR.pi() - z
     840            if w>0:
     841                z = z + RR.pi()
     842            z /= a
     843            return z
     844
     845        else:                    #Disconnected Case, disc > 0
     846            real_roots.sort() # increasing order
     847            real_roots.reverse() # decreasing order e1>e2>e3
     848            try:
     849                assert len(real_roots) == 3
     850            except:
     851                raise ValueError, ' no or more than one real root despite disc < 0'
     852            e1, e2, e3 = real_roots
     853            w1, w2 = E.period_lattice().basis(precision)
     854            a = (e1 - e3).sqrt()
     855            b = (e1 - e2).sqrt()
     856            on_egg = (x < e1)
     857
     858            # if P is on the "egg", replace it by P+T3
     859            # where T3=(e3,y3) is a 2-torsion point on
     860            # the egg coming from w2/2 on the lattice
     861
     862            if on_egg:
     863                y3 = -(a1*e3+a3)/2
     864                lam = (y-y3)/(x-e3)
     865                x3 = lam*(lam+a1)-a2-x-e3
     866                y = lam*(x-x3)-y-a1*x3-a3
     867                x = x3
     868            c = (x - e3).sqrt()
     869            while (a-b)/a > 0.5**precision:
     870                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     871
     872            z = (a/c).arcsin()/a
     873            if (2*y+a1*x+a3) > 0:
     874                z = w1 - z
     875            if on_egg:
     876                z = z + w2/2
     877            return z
     878
     879    ##############################  end  ################################
     880
    732881
    733882class EllipticCurvePoint_finite_field(EllipticCurvePoint_field):
    734883    def _magma_init_(self):
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r 91af4c3f6b92 -r 60d3d284fd13 sage/schemes/elliptic_curves/ell_rational_field.py
    a b AUTHORS: 
    1111   -- Christian Wuthrich (2007): added padic sha computation
    1212   -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files.
    1313   -- John Cremona (2008-01)
     14   -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points
     15   -- John Cremona (2008-07): further work on integral_points
    1416"""
    1517
    1618#*****************************************************************************
    import sage.databases.cremona 
    5557import sage.databases.cremona
    5658from   sage.libs.pari.all import pari
    5759import sage.functions.transcendental as transcendental
    58 import math
     60from sage.calculus.calculus import sqrt, floor, ceil
    5961import sage.libs.mwrank.all as mwrank
    6062import constructor
    6163from sage.interfaces.all import gp
    import ell_tate_curve 
    8486import ell_tate_curve
    8587
    8688factor = arith.factor
    87 sqrt = math.sqrt
    88 exp = math.exp
    8989mul = misc.mul
    9090next_prime = arith.next_prime
    9191kronecker_symbol = arith.kronecker_symbol
    Q = RationalField() 
    9393Q = RationalField()         
    9494C = ComplexField()
    9595R = RealField()
     96Z = IntegerRing()
    9697IR = rings.RealIntervalField(20)
    9798
    9899_MAX_HEIGHT=21
    class EllipticCurve_rational_field(Ellip 
    14901491        self.__regulator[proof] = reg
    14911492        return self.__regulator[proof]
    14921493
     1494    def height_pairing_matrix(self, points=None, precision=None):
     1495        """
     1496        Returns the height pairing matrix of the given points on this
     1497        curve, which must be defined over Q.
     1498
     1499        INPUT:
     1500
     1501        points -- either a list of points, which must be on this
     1502                  curve, or (default) None, in which case self.gens()
     1503                  will be used.
     1504        precision -- number of bit of precision of result
     1505                 (default: default RealField)
     1506                 
     1507        TODO: implement variable precision for heights of points
     1508
     1509        EXAMPLES:
     1510            sage: E = EllipticCurve([0, 0, 1, -1, 0])
     1511            sage: E.height_pairing_matrix()
     1512            [0.0511114082399688]
     1513           
     1514            For rank 0 curves, the result is a valid 0x0 matrix:
     1515            sage: EllipticCurve('11a').height_pairing_matrix()
     1516            []
     1517            sage: E=EllipticCurve('5077a1')
     1518            sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]])             
     1519            [  1.36857250535393  -1.30957670708658 -0.634867157837156]
     1520            [ -1.30957670708658   2.71735939281229   1.09981843056673]
     1521            [-0.634867157837156   1.09981843056673  0.668205165651928]
     1522
     1523        """
     1524        if points is None:
     1525            points = self.gens()
     1526        else:
     1527            for P in points:
     1528                assert P.curve() is self
     1529
     1530        r = len(points)
     1531        if precision is None:
     1532            RR = rings.RealField()
     1533        else:
     1534            RR = rings.RealField(precision)
     1535        M = matrix.MatrixSpace(RR, r)
     1536        mat = M()
     1537        for j in range(r):
     1538            mat[j,j] = points[j].height()
     1539        for j in range(r):
     1540            for k in range(j+1,r):             
     1541                mat[j,k]=((points[j]+points[k]).height() - mat[j,j] - mat[k,k])/2
     1542                mat[k,j]=mat[j,k]
     1543        return mat
     1544
     1545
    14931546    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False):
    14941547        """
    14951548        Given a list of rational points on E, compute the saturation
    class EllipticCurve_rational_field(Ellip 
    23102363        OUTPUT:
    23112364            The EllipticCurveTorsionSubgroup instance associated to this elliptic curve.
    23122365
     2366        NOTE:
     2367            To see the torsion points as a list, use torsion_points()
     2368
    23132369        EXAMPLES:
    23142370            sage: EllipticCurve('11a').torsion_subgroup()
    23152371            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 
    23362392            self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag)
    23372393            self.__torsion_order = self.__torsion_subgroup.order()
    23382394            return self.__torsion_subgroup
     2395
     2396    def torsion_points(self, flag=0):
     2397        """
     2398        Returns the torsion points of this elliptic curve as a sorted list.
     2399
     2400        INPUT:
     2401            flag -- (default: 0)  chooses PARI algorithm:
     2402              flag = 0: uses Doud algorithm
     2403              flag = 1: uses Lutz-Nagell algorithm
     2404
     2405        OUTPUT:
     2406            A list of all the torsion points on this elliptic curve.
     2407
     2408        EXAMPLES:
     2409            sage: EllipticCurve('11a').torsion_points()
     2410            [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)]
     2411            sage: EllipticCurve('37b').torsion_points() 
     2412            [(0 : 1 : 0), (8 : -19 : 1), (8 : 18 : 1)]
     2413
     2414            sage: E=EllipticCurve([-1386747,368636886])   
     2415            sage: E.torsion_subgroup()
     2416            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
     2417            sage: E.torsion_points() 
     2418            [(-1293 : 0 : 1),
     2419            (-933 : -29160 : 1),
     2420            (-933 : 29160 : 1),
     2421            (-285 : -27216 : 1),
     2422            (-285 : 27216 : 1),
     2423            (0 : 1 : 0),
     2424            (147 : -12960 : 1),
     2425            (147 : 12960 : 1),
     2426            (282 : 0 : 1),
     2427            (1011 : 0 : 1),
     2428            (1227 : -22680 : 1),
     2429            (1227 : 22680 : 1),
     2430            (2307 : -97200 : 1),
     2431            (2307 : 97200 : 1),
     2432            (8787 : -816480 : 1),
     2433            (8787 : 816480 : 1)]
     2434        """
     2435        multiples = sage.groups.generic.multiples
     2436        gens = self.torsion_subgroup(flag).gens()
     2437        ngens = len(gens)
     2438        if ngens == 0:
     2439            return [self(0)]
     2440        pts = list(multiples(gens[0],gens[0].order()))
     2441        if ngens == 1:
     2442            pts.sort()
     2443            return pts
     2444        pts = [P+Q for P in pts for Q in multiples(gens[1],gens[1].order())]
     2445        pts.sort()
     2446        return pts
    23392447
    23402448    ## def newform_eval(self, z, prec):
    23412449##         """
    class EllipticCurve_rational_field(Ellip 
    35493657        Eq = ell_tate_curve.TateCurve(self,p)
    35503658        self._tate_curve[p] = Eq
    35513659        return Eq
     3660
     3661    def height(self, precision=53):
     3662        """Returns real height of this elliptic curve
     3663        This is used in integral_points()
     3664
     3665
     3666        INPUT: precision (integer) bit precision of the field of real
     3667               numbers in which the result will lie (default: 53 as in
     3668               RealField())
     3669       
     3670        EXAMPLES:
     3671            sage: E=EllipticCurve('5077a1')
     3672            sage: E.height()               
     3673            17.4513334798896
     3674            sage: E.height(100)
     3675            17.451333479889612702508579399
     3676            sage: E=EllipticCurve([0,0,0,0,1])
     3677            sage: E.height()                 
     3678            1.38629436111989
     3679            sage: E=EllipticCurve([0,0,0,1,0])
     3680            sage: E.height()                 
     3681            7.45471994936400
     3682
     3683        """
     3684        #Computes height of curve 'self' according to Co1
     3685        R = RealField(precision)
     3686        c4 = self.c4()
     3687        c6 = self.c6()
     3688        j = self.j_invariant()
     3689        log_g2 = R((c4/12)).abs().log()
     3690        log_g3 = R((c6/216)).abs().log()
     3691
     3692        if j == 0:
     3693            h_j = R(1)
     3694        else:
     3695            h_j = max(log(R(abs(j.numerator()))), log(R(abs(j.denominator()))))
     3696        if (self.c4() != 0) and (self.c6() != 0):
     3697            h_gs = max(1, log_g2, log_g3)
     3698        elif c4 == 0:
     3699            if c6 == 0:
     3700                return max(1,h_j)
     3701            h_gs = max(1, log_g3)
     3702        else:
     3703            h_gs = max(1, log_g2)
     3704        return max(R(1),h_j, h_gs)
     3705
     3706    def integral_points(self, mw_base='auto', both_signs=False):
     3707        """
     3708        Computes all integral points (up to sign) on the elliptic
     3709        curve E which has Mordell-Weil basis mw_base.
     3710       
     3711        INPUT:
     3712            self -- EllipticCurve_Rational_Field
     3713            mw_base -- list of EllipticCurvePoint generating the
     3714                       Mordell-Weil group of E
     3715                    (default: 'auto' - calls self.gens())
     3716
     3717            both_signs -- True/False (default False): if True the
     3718                       output contains both P and -P, otherwise only
     3719                       one of ecah pair.
     3720           
     3721        OUTPUT:
     3722            set of all integral points on E
     3723            (up to sign unless both_signs is True)
     3724           
     3725        HINTS:
     3726            - The complexity increases exponentially in the rank of curve E.
     3727            - It can help if you try another Mordell-Weil base, because the
     3728            computation time depends on this, too.   
     3729           
     3730        EXAMPLES:
     3731        A curve of rank 3 with no torsion points
     3732       
     3733            sage: E=EllipticCurve([0,0,1,-7,6])
     3734            sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
     3735            sage: a=E.integral_points([P1,P2,P3]); a 
     3736            [(-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)]
     3737           
     3738        It is also possible to not specify mw_base and tors_base but
     3739        then the Mordell-Weil basis must be found by other functions
     3740        which will takes longer
     3741       
     3742            sage: E=EllipticCurve([0,0,1,-7,6])
     3743            sage: a=E.integral_points(both_signs=True); a
     3744            [(-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)]
     3745           
     3746
     3747        An example with negative discriminant:
     3748            sage: EllipticCurve('900d1').integral_points()
     3749            [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)]
     3750       
     3751        Another example with rank 5 and no torsion points
     3752           
     3753            sage: E=EllipticCurve([-879984,319138704])
     3754            sage: P1=E.point((540,1188)); P2=E.point((576,1836))
     3755            sage: P3=E.point((468,3132)); P4=E.point((612,3132))
     3756            sage: P5=E.point((432,4428))
     3757            sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (400s!)
     3758            54
     3759           
     3760           
     3761        NOTES:
     3762            - This function uses the algorithm given in [Co1]
     3763            REFERENCES:
     3764                -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations
     3765                        GTM 239, Springer 2007
     3766        AUTHORS:
     3767            - Michael Mardaus (2008-07)
     3768            - Tobias Nagel (2008-07)
     3769            - John Cremona (2008-07)
     3770        """
     3771        #####################################################################
     3772        # INPUT CHECK #######################################################
     3773        try:
     3774            assert self.is_integral()
     3775        except:
     3776            raise ValueError, "integral_points() can only be called on an integral model"
     3777       
     3778        if mw_base=='auto':
     3779            mw_base = self.gens()
     3780            r = len(mw_base)
     3781        else:
     3782            try:
     3783                r = len(mw_base)
     3784            except TypeError:
     3785                raise TypeError, 'mw_base must be a list'
     3786            for P in mw_base:
     3787                try:
     3788                    assert P.curve() is self
     3789                except:
     3790                    raise ValueError, "points are not on the correct curve"
     3791               
     3792        tors_base = self.torsion_subgroup().gens()
     3793        len_tors = len(tors_base)
     3794        tors_points = self.torsion_points()
     3795
     3796        # END INPUT-CHECK####################################################
     3797        #####################################################################
     3798               
     3799        #####################################################################
     3800        # INTERNAL FUNCTIONS ################################################
     3801   
     3802        ############################## begin ################################
     3803        def point_preprocessing(list):
     3804            #Transforms mw_base so that at most one point is on the
     3805            #compact component of the curve
     3806            Q = []
     3807            mem = -1
     3808            for i in range(0,len(list)):
     3809                if not list[i].is_on_identity_component(): # i.e. is on "egg"
     3810                    if mem == -1:
     3811                        mem = i
     3812                    else:
     3813                        Q.append(list[i]+list[mem])
     3814                        mem = i
     3815                else:
     3816                    Q.append(list[i])
     3817            if mem != -1: #add last point, which is not in egg, to Q
     3818                Q.append(2*list[mem])
     3819            return Q
     3820        ##############################  end  ################################
     3821        ############################## begin ################################
     3822        def modified_height(i):#computes modified height if base point i
     3823            return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2)
     3824        ##############################  end  ################################
     3825        ############################## begin ################################
     3826        def search_points(xmin,xmax):
     3827            """Returns a list of integral points on self with x in the interval
     3828            """
     3829            if both_signs:
     3830                return sum([self.lift_x(x,all=True) for x in range(xmin,xmax)],[])
     3831            else:
     3832                return sum([self.lift_x(x,all=True)[1:] for x in range(xmin,xmax)],[])
     3833        ##############################  end  ################################
     3834        ############################## begin ################################
     3835        def search_remaining_points():
     3836            """Returns list of integral points on curve E written as
     3837               linear combination of n times the mordell-weil base
     3838               points and torsion points (n is bounded by H_q, which
     3839               will be computed at runtime)
     3840            """
     3841            OP=self(0)
     3842            pts=[]
     3843            for i in range(1,((2*H_q+1)**r)/2):
     3844                koeffs = Z(i).digits(base=2*H_q+1)
     3845                koeffs = [0]*(r-len(koeffs)) + koeffs
     3846                P = self(0)
     3847                for index in range(r):
     3848                    P = P + (koeffs[index]-H_q)*mw_base[index]
     3849                for Q in tors_points: # P + torsion points (includes 0)
     3850                    P = P + Q
     3851                    if P!=OP and P[0].is_integral():
     3852                        pts.append(P)
     3853                        if both_signs:
     3854                            Q = -P
     3855                            if not Q in pts:
     3856                                pts.append(Q)
     3857            return pts
     3858
     3859        ##############################  end  #################################
     3860   
     3861        # END Internal functions #############################################
     3862        ######################################################################
     3863   
     3864        if (r == 0):
     3865            pts = self.torsion_points()
     3866            pts = [P for P in pts if P[0].is_integral()]
     3867            pts.sort()
     3868            return pts
     3869
     3870        e = R(1).exp()
     3871        pi = R(constants.pi)
     3872   
     3873        g2 = self.c4()/12
     3874        g3 = self.c6()/216
     3875        disc = self.discriminant()
     3876        j = self.j_invariant()
     3877        b2 = self.b2()
     3878       
     3879        Qx = rings.PolynomialRing(RationalField(),'x')
     3880        pol = Qx([-self.c6()/216,-self.c4()/12,0,4])
     3881        if disc > 0: # two real component -> 3 roots in RR
     3882            e1,e2,e3 = pol.roots(R,multiplicities=False)
     3883            mw_base = point_preprocessing(mw_base) #at most one point in E^{egg}
     3884   
     3885        elif disc < 0: # one real component => 1 root in RR (=: e3),
     3886                       # 2 roots in C (e1,e2)
     3887            roots = pol.roots(C,multiplicities=False)
     3888            e3 = pol.roots(R,multiplicities=False)[0]
     3889            roots.remove(e3)
     3890            e1,e2 = roots
     3891               
     3892        # Algorithm presented in [Co1]
     3893        h_E = self.height()
     3894        w1, w2 = self.period_lattice().basis()
     3895        mu = R(disc).abs().log()
     3896        if j!=0:
     3897            mu += max(R(1),R(j).abs().log())
     3898        if b2!=0:
     3899            mu += max(R(1),R(b2).abs().log())
     3900            mu += log(R(2))
     3901       
     3902        c1 = (mu + 2.14).exp()
     3903        M = self.height_pairing_matrix(mw_base)
     3904        c2 = min(M.charpoly ().roots(multiplicities=False))
     3905        c3 = (w1**2)*R(b2).abs()/48 + 8
     3906        c5 = (c1*c3).sqrt()
     3907        c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
     3908       
     3909        mw_base_log = [] #contains \Phi(Q_i)
     3910        mod_h_list = []  #contains h_m(Q_i)
     3911        c9_help_list = []
     3912        for i in range(0,r):
     3913            mw_base_log.append(mw_base[i].elliptic_logarithm().abs())
     3914            mod_h_list.append(modified_height(i))
     3915            c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i])
     3916   
     3917        c8 = max(e*h_E,max(mod_h_list))
     3918        c9 = e/c7.sqrt() * min(c9_help_list)
     3919        n=r+1
     3920        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))
     3921   
     3922        top = 128 #arbitrary first upper bound
     3923        bottom = 0
     3924        log_c9=log(c9); log_c5=log(c5)
     3925        log_r_top = log(R(r*(10**top)))
     3926   
     3927        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):
     3928            #initial bound 'top' too small, upshift of search interval
     3929            bottom = top
     3930            top = 2*top
     3931        while top >= bottom: #binary-search like search for fitting exponent (bound)
     3932            bound = floor(bottom + (top - bottom)/2)
     3933            log_r_bound = log(R(r*(10**bound)))
     3934            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):
     3935                bottom = bound + 1
     3936            else:
     3937                top = bound - 1
     3938   
     3939        H_q = R(10)**bound
     3940        break_cond = 0 #at least one reduction step
     3941        #reduction via LLL
     3942        while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10%
     3943            c = R((H_q**n)*10)  #c has to be greater than H_q^n
     3944            M = matrix.MatrixSpace(Z,n)
     3945            m = M.identity_matrix()
     3946            for i in range(r):
     3947                m[i, r] = R(c*mw_base_log[i]).round()
     3948            m[r,r] = R(c*w1).round()
     3949   
     3950            #LLL - implemented in sage - operates on rows not on columns
     3951            m_LLL = m.LLL()
     3952            m_gram = m_LLL.gram_schmidt()[0]
     3953            b1_norm = R(m_LLL.row(0).norm())
     3954   
     3955            #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0
     3956            c1_LLL = -1
     3957            for i in range(n):
     3958                tmp = R(b1_norm/(m_gram.row(i).norm()))
     3959                if tmp > c1_LLL:
     3960                    c1_LLL = tmp
     3961   
     3962            if c1_LLL < 0:
     3963                raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base'
     3964           
     3965            d_L_0 = R(b1_norm**2 / c1_LLL)
     3966   
     3967            #Reducing of upper bound
     3968            Q = r * H_q**2
     3969            T = (1 + (3/2*r*H_q))/2
     3970            if d_L_0 < R(T**2+Q):
     3971                d_L_0 = 10*(T**2*Q)
     3972            low_bound = R((sqrt(d_L_0 - Q) - T)/c)
     3973   
     3974            #new bound according to low_bound and upper bound [c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
     3975            H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2)))
     3976            H_q_new = ceil(H_q_new)
     3977            break_cond = R(H_q_new/H_q)
     3978            H_q = H_q_new
     3979            #END LLL-Reduction loop
     3980   
     3981   
     3982        int_points = []
     3983        ##for a more detailed view how the function works uncomment
     3984        ##all print - statements below
     3985        if disc > 0:
     3986            ##Points in egg have e2>=x>=e1
     3987            int_points += search_points(ceil(e1), floor(e2)+1)
     3988            #print 'Points on compact component \n',list(int_points)
     3989           
     3990        ##Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3
     3991        int_points += search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3))))
     3992        #print 'Points on compact component and non-compact component part 1 \n', list(int_points)
     3993       
     3994        #print 'starting search of remainig points using bound ',H_q
     3995        int_points += search_remaining_points()
     3996        #print 'Integral points:\n',list(int_points)
     3997        #print 'Total amount of points:',len(int_points)
     3998        if both_signs:
     3999            int_points = list(set(int_points)) # remove duplicates
     4000        else:
     4001            xlist = set([P[0] for P in int_points])
     4002            int_points = [self.lift_x(x) for x in xlist]
     4003        int_points.sort()
     4004        return int_points
    35524005
    35534006
    35544007def cremona_curves(conductors):