Ticket #3674: sage-trac3674a.patch

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

Based on 3.0.4

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1216537839 -3600
    # Node ID ada4f959e6e5df25190883d9220ed7828be7f2f8
    # Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
    #3674: preliminary version
    
    diff -r 91af4c3f6b92 -r ada4f959e6e5 sage/schemes/elliptic_curves/ell_point.py
    a b class EllipticCurvePoint_field(AdditiveG 
    708708        h = self.curve().pari_curve().ellheight([self[0], self[1]])
    709709        return rings.RR(h)
    710710   
     711    def is_on_identity_component(self):
     712        """
     713        Returns True iff this point is on the identity component of
     714        its curve (over an ordered field)
     715
     716        INPUT:
     717            self -- a point on a curve over any ordered field (e.g. QQ)
     718
     719        OUTPUT:
     720            True iff the point is on the identity component of the identity
     721
     722        EXAMPLES:
     723            sage: E=EllipticCurve('5077a1')     
     724            sage: [E.lift_x(x).is_on_identity_component() for x in range(-3,5)]
     725            [False, False, False, False, False, True, True, True]
     726        """
     727        if self.is_zero():
     728            return True
     729        E = self.curve()
     730        if E.discriminant() < 0:
     731            return False
     732        gx =E.division_polynomial(2)
     733        gxd = gx.derivative()
     734        gxdd = gxd.derivative()
     735        return ( gxd(self[0]) > 0 and gxdd(self[0]) > 0)
     736
    711737    def xy(self):
    712738        """
    713739        Return the x and y coordinates of this point, as a 2-tuple.
    class EllipticCurvePoint_field(AdditiveG 
    729755            return self[0], self[1]
    730756        else:
    731757            return self[0]/self[2], self[1]/self[2]
     758
     759    def elliptic_logarithm(self, precision=53):
     760        """ Returns the elliptic logarithm of this point on an
     761            elliptic curve defined over the reals
     762
     763        INPUT: - precision: a positive integer (default 53) setting
     764        the number of bits of precision required
     765
     766        ALGORITHM: See -[Co2] Cohen H., A Course in Computational
     767                        Algebraic Number Theory GTM 138, Springer 1996
     768
     769        AUTHORS:
     770            - Michael Mardaus (2008-07) }
     771            - Tobias Nagel (2008-07)    } original version from [Co2]
     772            - John Cremona (2008-07)    revision following eclib code
     773       
     774        """
     775        RR = rings.RealField(precision)
     776        CC = rings.ComplexField(precision)
     777        if self.is_zero():
     778            return CC(0)
     779
     780        #Initialize
     781
     782        E = self.curve()
     783        a1 = RR(E.a1())
     784        a2 = RR(E.a2())
     785        a3 = RR(E.a3())
     786        b2 = RR(E.b2())
     787        b4 = RR(E.b4())
     788        disc = E.discriminant()
     789        x = RR(self[0])
     790        y = RR(self[1])
     791        real_roots = E.division_polynomial(2).roots(RR,multiplicities=False)
     792
     793        if disc < 0: #Connected Case
     794            # Here we use formulae equivalent to those in Cohen, but better
     795            # behaved when roots are close together
     796            try:
     797                assert len(real_roots) == 1
     798            except:
     799                raise ValueError, ' no or more than one real root despite disc < 0'
     800            e1 = real_roots[0]
     801            zz = (e1-e2).sqrt() # complex
     802            beta = (e1-e2).abs()
     803            a = 2*zz.abs()
     804            b = 2*zz.real();
     805            c = (x-e1+beta)/((x-e1).sqrt())
     806            while (a - b)/a > 0.5**precision:
     807                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     808
     809            z = (a/c).arcsin()
     810            w = 2*y+a1*x+a3
     811            if w*((x-e1)*(x-e1)-beta*beta) >= 0:
     812                z = pi - z
     813            if w>0:
     814                z = z + pi
     815            z /= a
     816            return z
     817
     818        else:                    #Disconnected Case, disc > 0
     819            real_roots.sort() # increasing order
     820            real_roots.reverse() # decreasing order e1>e2>e3
     821            try:
     822                assert len(real_roots) == 3
     823            except:
     824                raise ValueError, ' no or more than one real root despite disc < 0'
     825            e1, e2, e3 = real_roots
     826            w1, w2 = E.period_lattice().basis(precision)
     827            a = (e1 - e3).sqrt()
     828            b = (e1 - e2).sqrt()
     829            on_egg = (x < e1)
     830
     831            # if P is on the "egg", replace it by P+T3
     832            # where T3=(e3,y3) is a 2-torsion point on
     833            # the egg coming from w2/2 on the lattice
     834
     835            if on_egg:
     836                y3 = -(a1*e3+a3)/2
     837                lam = (y-y3)/(x-e3)
     838                x3 = lam*(lam+a1)-a2-x-e3
     839                y = lam*(x-x3)-y-a1*x3-a3
     840                x = x3
     841            c = (x - e3).sqrt()
     842            while (a-b)/a > 0.5**precision:
     843                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     844
     845            z = (a/c).arcsin()/a
     846            if (2*y+a1*x+a3) > 0:
     847                z = w1 - z
     848            if on_egg:
     849                z = z + w2/2
     850            return z
     851
     852    ##############################  end  ################################
     853
    732854
    733855class EllipticCurvePoint_finite_field(EllipticCurvePoint_field):
    734856    def _magma_init_(self):
  • sage/schemes/elliptic_curves/ell_rational_field.py

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