Ticket #3674: sage-trac3674new.patch

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

Replaces ALL above patches

  • sage/schemes/elliptic_curves/ell_generic.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1218142976 -3600
    # Node ID 00987cd79c91861457154b22b1516357b55cb158
    # Parent  5e2426d277ebfd4f7010fb71fee7d5ab3f3574b2
    #3674: patch implementing integral points over Q and related functions
    
    diff -r 5e2426d277eb -r 00987cd79c91 sage/schemes/elliptic_curves/ell_generic.py
    a b class EllipticCurve_generic(plane_curve. 
    497497            args = tuple(args[0])
    498498        return plane_curve.ProjectiveCurve_generic.__call__(self, *args, **kwds)
    499499   
     500    def is_x_coord(self, x):
     501        """
     502        Returns whether x is the x-coordinate of a point on this curve.
     503
     504        See also lift_x() to find the point(s) with a given
     505        x_coordinate.  This function may be useful in cases where
     506        testing an element of tha base field for being a square is
     507        faster than finding its square root.
     508       
     509        EXAMPLES:
     510            sage: E = EllipticCurve('37a'); E
     511            Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
     512            sage: E.is_x_coord(1)
     513            True
     514            sage: E.is_x_coord(2)
     515            True
     516           
     517        There are no rational points with x-cordinate 3.
     518            sage: E.is_x_coord(3)
     519            False
     520
     521        However, there are such points in $E(\R)$:
     522            sage: E.change_ring(RR).is_x_coord(3)
     523            True
     524           
     525        And of course it always works in $E(\C)$:
     526            sage: E.change_ring(RR).is_x_coord(-3)
     527            False
     528            sage: E.change_ring(CC).is_x_coord(-3)
     529            True
     530                                 
     531        AUTHOR: John Cremona, 2008-08-07 (adapted from lift_x())
     532           
     533        TEST:
     534            sage: E=EllipticCurve('5077a1')                     
     535            sage: [x for x in srange(-10,10) if E.is_x_coord (x)]
     536            [-3, -2, -1, 0, 1, 2, 3, 4, 8]
     537
     538            sage: F=GF(32,'a')
     539            sage: E=EllipticCurve(F,[1,0,0,0,1])
     540            sage: set([P[0] for P in E.points() if P!=E(0)]) == set([x for x in F if E.is_x_coord(x)])
     541            True
     542
     543        """
     544        K = self.base_ring()
     545        try:
     546            x = K(x)
     547        except TypeError:
     548            raise TypeError, 'x must be coercible into the base ring of the curve'
     549        a1, a2, a3, a4, a6 = self.ainvs()
     550        fx = ((x + a2) * x + a4) * x + a6
     551        if a1.is_zero() and a3.is_zero():
     552            return fx.is_square()
     553        b = (a1*x + a3)
     554        if K.characteristic() == 2:
     555            R = PolynomialRing(K, 'y')
     556            F = R([-fx,b,1])
     557            return len(F.roots())>0
     558        D = b*b + 4*fx
     559        return D.is_square()
     560
    500561    def lift_x(self, x, all=False):
    501562        """
    502563        Given the x-coordinate of a point on the curve, use the defining
  • sage/schemes/elliptic_curves/ell_point.py

    diff -r 5e2426d277eb -r 00987cd79c91 sage/schemes/elliptic_curves/ell_point.py
    a 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-1):
     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 5e2426d277eb -r 00987cd79c91 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 
    22982351            self.__torsion_order = self.torsion_subgroup().order()
    22992352            return self.__torsion_order
    23002353
    2301     def torsion_subgroup(self, flag=0):
     2354    def torsion_subgroup(self, algorithm="pari"):
    23022355        """
    23032356        Returns the torsion subgroup of this elliptic curve.
    23042357
    23052358        INPUT:
    2306             flag -- (default: 0)  chooses PARI algorithm:
    2307               flag = 0: uses Doud algorithm
    2308               flag = 1: uses Lutz-Nagell algorithm
     2359            algorithm -- string:
     2360                 "pari"  -- (default) use the pari library
     2361                 "doud" -- use Doud's algorithm
     2362                 "lutz_nagell" -- use the Lutz-Nagell theorem
    23092363
    23102364        OUTPUT:
    23112365            The EllipticCurveTorsionSubgroup instance associated to this elliptic curve.
     2366
     2367        NOTE:
     2368            To see the torsion points as a list, use torsion_points()
    23122369
    23132370        EXAMPLES:
    23142371            sage: EllipticCurve('11a').torsion_subgroup()
    class EllipticCurve_rational_field(Ellip 
    23332390        try:
    23342391            return self.__torsion_subgroup
    23352392        except AttributeError:
    2336             self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag)
     2393            self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, algorithm)
    23372394            self.__torsion_order = self.__torsion_subgroup.order()
    23382395            return self.__torsion_subgroup
     2396
     2397    def torsion_points(self, algorithm="pari"):
     2398        """
     2399        Returns the torsion points of this elliptic curve as a sorted list.
     2400
     2401        INPUT:
     2402            algorithm -- string:
     2403                 "pari"  -- (default) use the pari library
     2404                 "doud" -- use Doud's algorithm
     2405                 "lutz_nagell" -- use the Lutz-Nagell theorem
     2406
     2407        OUTPUT:
     2408            A list of all the torsion points on this elliptic curve.
     2409
     2410        EXAMPLES:
     2411            sage: EllipticCurve('11a').torsion_points()
     2412            [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)]
     2413            sage: EllipticCurve('37b').torsion_points() 
     2414            [(0 : 1 : 0), (8 : -19 : 1), (8 : 18 : 1)]
     2415
     2416            sage: E=EllipticCurve([-1386747,368636886])   
     2417            sage: T=E.torsion_subgroup(); T
     2418            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
     2419            sage: T == E.torsion_subgroup(algorithm="doud")
     2420            True
     2421            sage: T == E.torsion_subgroup(algorithm="lutz_nagell")
     2422            True
     2423            sage: E.torsion_points() 
     2424            [(-1293 : 0 : 1),
     2425            (-933 : -29160 : 1),
     2426            (-933 : 29160 : 1),
     2427            (-285 : -27216 : 1),
     2428            (-285 : 27216 : 1),
     2429            (0 : 1 : 0),
     2430            (147 : -12960 : 1),
     2431            (147 : 12960 : 1),
     2432            (282 : 0 : 1),
     2433            (1011 : 0 : 1),
     2434            (1227 : -22680 : 1),
     2435            (1227 : 22680 : 1),
     2436            (2307 : -97200 : 1),
     2437            (2307 : 97200 : 1),
     2438            (8787 : -816480 : 1),
     2439            (8787 : 816480 : 1)]
     2440        """
     2441        multiples = sage.groups.generic.multiples
     2442        gens = self.torsion_subgroup(algorithm).gens()
     2443        ngens = len(gens)
     2444        if ngens == 0:
     2445            return [self(0)]
     2446        pts = list(multiples(gens[0],gens[0].order()))
     2447        if ngens == 1:
     2448            pts.sort()
     2449            return pts
     2450        pts = [P+Q for P in pts for Q in multiples(gens[1],gens[1].order())]
     2451        pts.sort()
     2452        return pts
    23392453
    23402454    ## def newform_eval(self, z, prec):
    23412455##         """
    class EllipticCurve_rational_field(Ellip 
    33693483        we return v = -1.
    33703484
    33713485        INPUT:
    3372             D (int) -- (deault: 0) Heegner discriminant; if 0, use the
     3486            D (int) -- (default: 0) Heegner discriminant; if 0, use the
    33733487                       first discriminant < -4 that satisfies the Heegner
    33743488                       hypothesis
    33753489            verbose (bool) -- (default: True)
    class EllipticCurve_rational_field(Ellip 
    35493663        Eq = ell_tate_curve.TateCurve(self,p)
    35503664        self._tate_curve[p] = Eq
    35513665        return Eq
     3666
     3667    def height(self, precision=53):
     3668        """
     3669        Returns the real height of this elliptic curve.
     3670        This is used in integral_points()
     3671
     3672
     3673        INPUT:
     3674
     3675            precision -- (integer: default 53) bit precision of the
     3676               field of real numbers in which the result will lie
     3677       
     3678        EXAMPLES:
     3679            sage: E=EllipticCurve('5077a1')
     3680            sage: E.height()               
     3681            17.4513334798896
     3682            sage: E.height(100)
     3683            17.451333479889612702508579399
     3684            sage: E=EllipticCurve([0,0,0,0,1])
     3685            sage: E.height()                 
     3686            1.38629436111989
     3687            sage: E=EllipticCurve([0,0,0,1,0])
     3688            sage: E.height()                 
     3689            7.45471994936400
     3690
     3691        """
     3692        R = RealField(precision)
     3693        c4 = self.c4()
     3694        c6 = self.c6()
     3695        j = self.j_invariant()
     3696        log_g2 = R((c4/12)).abs().log()
     3697        log_g3 = R((c6/216)).abs().log()
     3698
     3699        if j == 0:
     3700            h_j = R(1)
     3701        else:
     3702            h_j = max(log(R(abs(j.numerator()))), log(R(abs(j.denominator()))))
     3703        if (self.c4() != 0) and (self.c6() != 0):
     3704            h_gs = max(1, log_g2, log_g3)
     3705        elif c4 == 0:
     3706            if c6 == 0:
     3707                return max(1,h_j)
     3708            h_gs = max(1, log_g3)
     3709        else:
     3710            h_gs = max(1, log_g2)
     3711        return max(R(1),h_j, h_gs)
     3712
     3713    def antilogarithm(self, z, precision=100):
     3714        """
     3715        Returns the rational point (if any) associated to this complex
     3716        number; the inverse of the elliptic logarithm function.
     3717
     3718        INPUT:
     3719        z -- a complex number representing an element of
     3720             CC/L where L is the period lattice of the elliptic
     3721             curve
     3722
     3723        precision -- an integer specifying the precision (in bits) which
     3724                will be used for the computation
     3725
     3726        OUTPUT: The rational point which is the image of z under
     3727        the Weierstrass parametrization, if it exists and can
     3728        be determined from z with default precision.
     3729
     3730        NOTE: This uses the function ellztopoint from the pari library
     3731
     3732        TODO: Extend the wrapping of ellztopoint() to allow passing of
     3733        the precision parameter.
     3734           
     3735        """
     3736       
     3737        E_pari = self.pari_curve(prec)
     3738        try:
     3739            coords = E_pari.ellztopoint(z)
     3740            if len(coords) == 1:
     3741                return self(0)
     3742            return self([CC(xy).real().simplest_rational() for xy in coords])
     3743        except TypeError:
     3744            raise ValueError, "No rational point computable from z"
     3745
     3746    def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
     3747        """
     3748        Computes all integral points (up to sign) on this elliptic curve.
     3749       
     3750        INPUT:
     3751            mw_base -- list of EllipticCurvePoint generating the
     3752                       Mordell-Weil group of E
     3753                    (default: 'auto' - calls self.gens())
     3754            both_signs -- True/False (default False): if True the
     3755                       output contains both P and -P, otherwise only
     3756                       one of each pair.         
     3757            verbose -- True/False (default False): if True, some
     3758                       details of the computation are output
     3759                                   
     3760        OUTPUT:
     3761            A sorted list of all the integral points on E
     3762            (up to sign unless both_signs is True)
     3763           
     3764        NOTES: The complexity increases exponentially in the rank of
     3765            curve E.  The computation time (but not the output!)
     3766            depends on the Mordell-Weil basis.  If mw_base is given
     3767            but it not a basis for the Mordell-Weil group (modulo
     3768            torsion), integral points which are not in the subgroup
     3769            generated by the given points will almost certainly not be
     3770            listed.
     3771           
     3772        EXAMPLES:
     3773        A curve of rank 3 with no torsion points
     3774       
     3775            sage: E=EllipticCurve([0,0,1,-7,6])
     3776            sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
     3777            sage: a=E.integral_points([P1,P2,P3]); a 
     3778            [(-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)]
     3779           
     3780            sage: a = E.integral_points([P1,P2,P3], verbose=True)
     3781            Using mw_basis  [(2 : 0 : 1), (4 : 6 : 1), (114/49 : -720/343 : 1)]
     3782            e1,e2,e3:  -3.01243037259331 1.06582054769620 1.94660982489710
     3783            Minimal eigenvalue of height pairing matrix:  0.472730555831538
     3784            x-coords of points on compact component with  -3 <=x<= 1
     3785            set([0, -1, -3, -2, 1])
     3786            x-coords of points on non-compact component with  2 <=x<= 6
     3787            set([2, 3, 4])
     3788            starting search of remaining points using coefficient bound  6
     3789            x-coords of extra integral points:
     3790            set([2, 3, 4, 37, 406, 8, 11, 14, 816, 52, 21, 342, 93])
     3791            Total number of integral points: 18
     3792
     3793        It is also possible to not specify mw_base, but then the
     3794        Mordell-Weil basis must be computed;  this may take much longer
     3795       
     3796            sage: E=EllipticCurve([0,0,1,-7,6])
     3797            sage: a=E.integral_points(both_signs=True); a
     3798            [(-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)]
     3799           
     3800
     3801        An example with negative discriminant:
     3802            sage: EllipticCurve('900d1').integral_points()
     3803            [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)]
     3804       
     3805        Another example with rank 5 and no torsion points
     3806           
     3807            sage: E=EllipticCurve([-879984,319138704])
     3808            sage: P1=E.point((540,1188)); P2=E.point((576,1836))
     3809            sage: P3=E.point((468,3132)); P4=E.point((612,3132))
     3810            sage: P5=E.point((432,4428))
     3811            sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (400s!)
     3812            54
     3813           
     3814           
     3815        NOTES:
     3816            - This function uses the algorithm given in [Co1]
     3817            REFERENCES:
     3818                -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations
     3819                        GTM 239, Springer 2007
     3820        AUTHORS:
     3821            - Michael Mardaus (2008-07)
     3822            - Tobias Nagel (2008-07)
     3823            - John Cremona (2008-07)
     3824        """
     3825        #####################################################################
     3826        # INPUT CHECK #######################################################
     3827        if not self.is_integral():
     3828            raise ValueError, "integral_points() can only be called on an integral model"
     3829       
     3830        if mw_base=='auto':
     3831            mw_base = self.gens()
     3832            r = len(mw_base)
     3833        else:
     3834            try:
     3835                r = len(mw_base)
     3836            except TypeError:
     3837                raise TypeError, 'mw_base must be a list'
     3838            if not all([P.curve() is self for P in mw_base]):
     3839                raise ValueError, "points are not on the correct curve"
     3840               
     3841        tors_points = self.torsion_points()
     3842
     3843        # END INPUT-CHECK####################################################
     3844        #####################################################################
     3845               
     3846        #####################################################################
     3847        # INTERNAL FUNCTIONS ################################################
     3848   
     3849        ############################## begin ################################
     3850        def point_preprocessing(list):
     3851            #Transforms mw_base so that at most one point is on the
     3852            #compact component of the curve
     3853            Q = []
     3854            mem = -1
     3855            for i in range(0,len(list)):
     3856                if not list[i].is_on_identity_component(): # i.e. is on "egg"
     3857                    if mem == -1:
     3858                        mem = i
     3859                    else:
     3860                        Q.append(list[i]+list[mem])
     3861                        mem = i
     3862                else:
     3863                    Q.append(list[i])
     3864            if mem != -1: #add last point, which is not in egg, to Q
     3865                Q.append(2*list[mem])
     3866            return Q
     3867        ##############################  end  ################################
     3868        ############################## begin ################################
     3869        def modified_height(i):#computes modified height if base point i
     3870            return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2)
     3871        ##############################  end  ################################
     3872        ############################## begin ################################
     3873        def integral_x_coords_in_interval(xmin,xmax):
     3874            """
     3875            Returns the set of integers x with xmin<=x<=xmax which are
     3876            x-coordinates of points on this curve.
     3877            """
     3878            return set([x for x  in range(xmin,xmax) if self.is_x_coord(x)])
     3879        ##############################  end  ################################
     3880        ############################## begin ################################
     3881        def integral_points_with_bounded_mw_ceoffs():
     3882            r"""
     3883            Returns the set of integers x which are x-coordinates of
     3884            points on the curve which are linear combinations of the
     3885            generators (basis and torsion points) with coefficients
     3886            bounded by $H_q$.  The bound $H_q$ will be computed at
     3887            runtime.
     3888            """
     3889            xs=set()
     3890            for i in range(ceil(((2*H_q+1)**r)/2)):
     3891                koeffs = Z(i).digits(base=2*H_q+1)
     3892                koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s
     3893                P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0))
     3894                for Q in tors_points: # P + torsion points (includes 0)
     3895                    tmp = P + Q
     3896                    if not tmp.is_zero():
     3897                        x = tmp[0]
     3898                        if x.is_integral():
     3899                            xs.add(x)
     3900            return xs
     3901        ##############################  end  #################################
     3902   
     3903        # END Internal functions #############################################
     3904        ######################################################################
     3905   
     3906        if (r == 0):
     3907            int_points = [P for P in tors_points if not P.is_zero()]
     3908            int_points = [P for P in int_points if P[0].is_integral()]
     3909            if not both_signs:
     3910                xlist = set([P[0] for P in int_points])
     3911                int_points = [self.lift_x(x) for x in xlist]
     3912            int_points.sort()
     3913            if verbose:
     3914                print 'Total number of integral points:',len(int_points)
     3915            return int_points
     3916
     3917
     3918        g2 = self.c4()/12
     3919        g3 = self.c6()/216
     3920        disc = self.discriminant()
     3921        j = self.j_invariant()
     3922        b2 = self.b2()
     3923       
     3924        Qx = rings.PolynomialRing(RationalField(),'x')
     3925        pol = Qx([-self.c6()/216,-self.c4()/12,0,4])
     3926        if disc > 0: # two real component -> 3 roots in RR
     3927            #on curve 897e4, only one root is found with default precision!
     3928            RR = R
     3929            prec = RR.precision()
     3930            ei = pol.roots(RR,multiplicities=False)
     3931            while len(ei)<3:
     3932                prec*=2
     3933                RR=RealField(prec)
     3934                ei = pol.roots(RR,multiplicities=False)
     3935            e1,e2,e3 = ei
     3936            if r >= 2: #preprocessing of mw_base only necessary if rank > 1
     3937                mw_base = point_preprocessing(mw_base) #at most one point in
     3938                                                       #E^{egg}, saved in P_egg
     3939
     3940        elif disc < 0: # one real component => 1 root in RR (=: e3),
     3941                       # 2 roots in C (e1,e2)
     3942            roots = pol.roots(C,multiplicities=False)
     3943            e3 = pol.roots(R,multiplicities=False)[0]
     3944            roots.remove(e3)
     3945            e1,e2 = roots
     3946 
     3947        e = R(1).exp()
     3948        pi = R(constants.pi)
     3949   
     3950        if verbose:
     3951            print "Using mw_basis ",mw_base
     3952            print "e1,e2,e3: ",e1,e2,e3
     3953           
     3954        # Algorithm presented in [Co1]
     3955        h_E = self.height()
     3956        w1, w2 = self.period_lattice().basis()
     3957        mu = R(disc).abs().log() / 6
     3958        if j!=0:
     3959            mu += max(R(1),R(j).abs().log()) / 6
     3960        if b2!=0:
     3961            mu += max(R(1),R(b2).abs().log())
     3962            mu += log(R(2))
     3963        else:
     3964            mu += 1
     3965       
     3966        c1 = (mu + 2.14).exp()
     3967        M = self.height_pairing_matrix(mw_base)
     3968        c2 = min(M.charpoly ().roots(multiplicities=False))
     3969        if verbose:
     3970            print "Minimal eigenvalue of height pairing matrix: ", c2
     3971
     3972        c3 = (w1**2)*R(b2).abs()/48 + 8
     3973        c5 = (c1*c3).sqrt()
     3974        c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
     3975       
     3976        mw_base_log = [] #contains \Phi(Q_i)
     3977        mod_h_list = []  #contains h_m(Q_i)
     3978        c9_help_list = []
     3979        for i in range(0,r):
     3980            mw_base_log.append(mw_base[i].elliptic_logarithm().abs())
     3981            mod_h_list.append(modified_height(i))
     3982            c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i])
     3983   
     3984        c8 = max(e*h_E,max(mod_h_list))
     3985        c9 = e/c7.sqrt() * min(c9_help_list)
     3986        n=r+1
     3987        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))
     3988   
     3989        top = 128 #arbitrary first upper bound
     3990        bottom = 0
     3991        log_c9=log(c9); log_c5=log(c5)
     3992        log_r_top = log(R(r*(10**top)))
     3993#        if verbose:
     3994#            print "[bottom,top] = ",[bottom,top]
     3995           
     3996        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):
     3997            #initial bound 'top' too small, upshift of search interval
     3998            bottom = top
     3999            top = 2*top
     4000        while top >= bottom: #binary-search like search for fitting exponent (bound)
     4001#            if verbose:
     4002#                print "[bottom,top] = ",[bottom,top]
     4003            bound = floor(bottom + (top - bottom)/2)
     4004            log_r_bound = log(R(r*(10**bound)))
     4005            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):
     4006                bottom = bound + 1
     4007            else:
     4008                top = bound - 1
     4009   
     4010        H_q = R(10)**bound
     4011        break_cond = 0 #at least one reduction step
     4012        #reduction via LLL
     4013        while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10%
     4014            c = R((H_q**n)*10)  #c has to be greater than H_q^n
     4015            M = matrix.MatrixSpace(Z,n)
     4016            m = M.identity_matrix()
     4017            for i in range(r):
     4018                m[i, r] = R(c*mw_base_log[i]).round()
     4019            m[r,r] = R(c*w1).round()
     4020   
     4021            #LLL - implemented in sage - operates on rows not on columns
     4022            m_LLL = m.LLL()
     4023            m_gram = m_LLL.gram_schmidt()[0]
     4024            b1_norm = R(m_LLL.row(0).norm())
     4025   
     4026            #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0
     4027            c1_LLL = -1
     4028            for i in range(n):
     4029                tmp = R(b1_norm/(m_gram.row(i).norm()))
     4030                if tmp > c1_LLL:
     4031                    c1_LLL = tmp
     4032   
     4033            if c1_LLL < 0:
     4034                raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base'
     4035           
     4036            d_L_0 = R(b1_norm**2 / c1_LLL)
     4037   
     4038            #Reducing of upper bound
     4039            Q = r * H_q**2
     4040            T = (1 + (3/2*r*H_q))/2
     4041            if d_L_0 < R(T**2+Q):
     4042                d_L_0 = 10*(T**2*Q)
     4043            low_bound = R((sqrt(d_L_0 - Q) - T)/c)
     4044   
     4045            #new bound according to low_bound and upper bound
     4046            #[c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
     4047            if low_bound != 0:
     4048                H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2)))
     4049                H_q_new = ceil(H_q_new)
     4050                if H_q_new == 1:
     4051                    break_cond = 1 # stops reduction
     4052                else:
     4053                    break_cond = R(H_q_new/H_q)
     4054                H_q = H_q_new
     4055            else:
     4056                break_cond = 1 # stops reduction, so using last H_q > 0
     4057            #END LLL-Reduction loop
     4058   
     4059        b2_12 = b2/12
     4060        if disc > 0:
     4061            ##Points in egg have X(P) between e1 and e2 [X(P)=x(P)+b2/12]:
     4062            x_int_points = integral_x_coords_in_interval(ceil(e1-b2_12), floor(e2-b2_12)+1)
     4063            if verbose:
     4064                print 'x-coords of points on compact component with ',ceil(e1-b2_12),'<=x<=',floor(e2-b2_12)
     4065                print x_int_points
     4066        else:
     4067            x_int_points = set()
     4068           
     4069        ##Points in noncompact component with X(P)< 2*max(|e1|,|e2|,|e3|) , espec. X(P)>=e3
     4070        x0 = ceil(e3-b2_12)
     4071        x1 = ceil(2*max(abs(e1),abs(e2),abs(e3)) - b2_12)
     4072        x_int_points2 = integral_x_coords_in_interval(x0, x1)
     4073        x_int_points = x_int_points.union(x_int_points2)
     4074        if verbose:
     4075            print 'x-coords of points on non-compact component with ',x0,'<=x<=',x1-1
     4076            print x_int_points2
     4077       
     4078        if verbose:
     4079            print 'starting search of remaining points using coefficient bound ',H_q
     4080        x_int_points3 = integral_points_with_bounded_mw_ceoffs()
     4081        x_int_points = x_int_points.union(x_int_points3)
     4082        if verbose:
     4083            print 'x-coords of extra integral points:'
     4084            print x_int_points3
     4085
     4086        # Now we have all the x-coordinates of integral points, and we
     4087        # construct the points, depending on the parameter both_signs:
     4088        if both_signs:
     4089            int_points = sum([self.lift_x(x,all=True) for x in x_int_points],[])
     4090        else:
     4091            int_points = [self.lift_x(x) for x in x_int_points]
     4092        int_points.sort()
     4093        if verbose:
     4094            print 'Total number of integral points:',len(int_points)
     4095        return int_points
    35524096
    35534097
    35544098def cremona_curves(conductors):