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

File 3674-jcremona-integral-points.patch, 40.1 KB (added by ncalexan, 9 years ago)
  • 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
    * * *
    #3674: extra changes
    * * *
    #3674: some bug fixes + more efficient version of integral_points_with_bounded_mw_ceoffs()
    
    diff -r 354bcfd4fd38 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 354bcfd4fd38 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
     822            # For some curves (e.g. 3314b3) the default precision is not enough!
     823            while len(real_roots)!=1:
     824                precision*=2
     825                RR=rings.RealField(precision)
     826                real_roots = pol.roots(RR,multiplicities=False)
     827            try:
     828                assert len(real_roots) == 1
     829            except:
     830                raise ValueError, ' no or more than one real root despite disc < 0'
     831            e1 = real_roots[0]
     832            roots = pol.roots(rings.ComplexField(precision),multiplicities=False)
     833            roots.remove(e1)
     834            e2,e3 = roots
     835            zz = (e1-e2).sqrt() # complex
     836            beta = (e1-e2).abs()
     837            a = 2*zz.abs()
     838            b = 2*zz.real();
     839            c = (x-e1+beta)/((x-e1).sqrt())
     840            while (a - b)/a > 0.5**(precision-1):
     841                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     842            z = (a/c).arcsin()
     843            w = 2*y+a1*x+a3           
     844            if w*((x-e1)*(x-e1)-beta*beta) >= 0:
     845                z = RR.pi() - z
     846            if w>0:
     847                z = z + RR.pi()
     848            z /= a
     849            return z
     850
     851        else:                    #Disconnected Case, disc > 0
     852            # For some curves (e.g. 2370i5) the default precision is not enough!
     853            while len(real_roots)!=3:
     854                precision*=2
     855                RR=rings.RealField(precision)
     856                real_roots = pol.roots(RR,multiplicities=False)
     857
     858            real_roots.sort() # increasing order
     859            real_roots.reverse() # decreasing order e1>e2>e3
     860            try:
     861                assert len(real_roots) == 3
     862            except:
     863                raise ValueError, ' no or more than one real root despite disc < 0'
     864            e1, e2, e3 = real_roots
     865            w1, w2 = E.period_lattice().basis(precision)
     866            a = (e1 - e3).sqrt()
     867            b = (e1 - e2).sqrt()
     868            on_egg = (x < e1)
     869
     870            # if P is on the "egg", replace it by P+T3
     871            # where T3=(e3,y3) is a 2-torsion point on
     872            # the egg coming from w2/2 on the lattice
     873
     874            if on_egg:
     875                y3 = -(a1*e3+a3)/2
     876                lam = (y-y3)/(x-e3)
     877                x3 = lam*(lam+a1)-a2-x-e3
     878                y = lam*(x-x3)-y-a1*x3-a3
     879                x = x3
     880            c = (x - e3).sqrt()
     881            while (a - b)/a > 0.5**(precision-1):
     882                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
     883
     884            z = (a/c).arcsin()/a
     885            if (2*y+a1*x+a3) > 0:
     886                z = w1 - z
     887            if on_egg:
     888                z = z + w2/2
     889            return z
     890
     891    ##############################  end  ################################
     892
    732893
    733894class EllipticCurvePoint_finite_field(EllipticCurvePoint_field):
    734895    def _magma_init_(self):
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r 354bcfd4fd38 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 
    15051506        self.__regulator[proof] = reg
    15061507        return self.__regulator[proof]
    15071508
     1509    def height_pairing_matrix(self, points=None, precision=None):
     1510        """
     1511        Returns the height pairing matrix of the given points on this
     1512        curve, which must be defined over Q.
     1513
     1514        INPUT:
     1515
     1516        points -- either a list of points, which must be on this
     1517                  curve, or (default) None, in which case self.gens()
     1518                  will be used.
     1519        precision -- number of bit of precision of result
     1520                 (default: default RealField)
     1521                 
     1522        TODO: implement variable precision for heights of points
     1523
     1524        EXAMPLES:
     1525            sage: E = EllipticCurve([0, 0, 1, -1, 0])
     1526            sage: E.height_pairing_matrix()
     1527            [0.0511114082399688]
     1528           
     1529            For rank 0 curves, the result is a valid 0x0 matrix:
     1530            sage: EllipticCurve('11a').height_pairing_matrix()
     1531            []
     1532            sage: E=EllipticCurve('5077a1')
     1533            sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]])             
     1534            [  1.36857250535393  -1.30957670708658 -0.634867157837156]
     1535            [ -1.30957670708658   2.71735939281229   1.09981843056673]
     1536            [-0.634867157837156   1.09981843056673  0.668205165651928]
     1537
     1538        """
     1539        if points is None:
     1540            points = self.gens()
     1541        else:
     1542            for P in points:
     1543                assert P.curve() == self
     1544
     1545        r = len(points)
     1546        if precision is None:
     1547            RR = rings.RealField()
     1548        else:
     1549            RR = rings.RealField(precision)
     1550        M = matrix.MatrixSpace(RR, r)
     1551        mat = M()
     1552        for j in range(r):
     1553            mat[j,j] = points[j].height()
     1554        for j in range(r):
     1555            for k in range(j+1,r):             
     1556                mat[j,k]=((points[j]+points[k]).height() - mat[j,j] - mat[k,k])/2
     1557                mat[k,j]=mat[j,k]
     1558        return mat
     1559
     1560
    15081561    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False):
    15091562        """
    15101563        Given a list of rational points on E, compute the saturation
    class EllipticCurve_rational_field(Ellip 
    23132366            self.__torsion_order = self.torsion_subgroup().order()
    23142367            return self.__torsion_order
    23152368
    2316     def torsion_subgroup(self, flag=0):
     2369    def torsion_subgroup(self, algorithm="pari"):
    23172370        """
    23182371        Returns the torsion subgroup of this elliptic curve.
    23192372
    23202373        INPUT:
    2321             flag -- (default: 0)  chooses PARI algorithm:
    2322               flag = 0: uses Doud algorithm
    2323               flag = 1: uses Lutz-Nagell algorithm
     2374            algorithm -- string:
     2375                 "pari"  -- (default) use the pari library
     2376                 "doud" -- use Doud's algorithm
     2377                 "lutz_nagell" -- use the Lutz-Nagell theorem
    23242378
    23252379        OUTPUT:
    23262380            The EllipticCurveTorsionSubgroup instance associated to this elliptic curve.
     2381
     2382        NOTE:
     2383            To see the torsion points as a list, use torsion_points()
    23272384
    23282385        EXAMPLES:
    23292386            sage: EllipticCurve('11a').torsion_subgroup()
    class EllipticCurve_rational_field(Ellip 
    23482405        try:
    23492406            return self.__torsion_subgroup
    23502407        except AttributeError:
    2351             self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag)
     2408            self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, algorithm)
    23522409            self.__torsion_order = self.__torsion_subgroup.order()
    23532410            return self.__torsion_subgroup
     2411
     2412    def torsion_points(self, algorithm="pari"):
     2413        """
     2414        Returns the torsion points of this elliptic curve as a sorted list.
     2415
     2416        INPUT:
     2417            algorithm -- string:
     2418                 "pari"  -- (default) use the pari library
     2419                 "doud" -- use Doud's algorithm
     2420                 "lutz_nagell" -- use the Lutz-Nagell theorem
     2421
     2422        OUTPUT:
     2423            A list of all the torsion points on this elliptic curve.
     2424
     2425        EXAMPLES:
     2426            sage: EllipticCurve('11a').torsion_points()
     2427            [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)]
     2428            sage: EllipticCurve('37b').torsion_points() 
     2429            [(0 : 1 : 0), (8 : -19 : 1), (8 : 18 : 1)]
     2430
     2431            sage: E=EllipticCurve([-1386747,368636886])   
     2432            sage: T=E.torsion_subgroup(); T
     2433            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
     2434            sage: T == E.torsion_subgroup(algorithm="doud")
     2435            True
     2436            sage: T == E.torsion_subgroup(algorithm="lutz_nagell")
     2437            True
     2438            sage: E.torsion_points() 
     2439            [(-1293 : 0 : 1),
     2440            (-933 : -29160 : 1),
     2441            (-933 : 29160 : 1),
     2442            (-285 : -27216 : 1),
     2443            (-285 : 27216 : 1),
     2444            (0 : 1 : 0),
     2445            (147 : -12960 : 1),
     2446            (147 : 12960 : 1),
     2447            (282 : 0 : 1),
     2448            (1011 : 0 : 1),
     2449            (1227 : -22680 : 1),
     2450            (1227 : 22680 : 1),
     2451            (2307 : -97200 : 1),
     2452            (2307 : 97200 : 1),
     2453            (8787 : -816480 : 1),
     2454            (8787 : 816480 : 1)]
     2455        """
     2456        multiples = sage.groups.generic.multiples
     2457        gens = self.torsion_subgroup(algorithm).gens()
     2458        ngens = len(gens)
     2459        if ngens == 0:
     2460            return [self(0)]
     2461        pts = list(multiples(gens[0],gens[0].order()))
     2462        if ngens == 1:
     2463            pts.sort()
     2464            return pts
     2465        pts = [P+Q for P in pts for Q in multiples(gens[1],gens[1].order())]
     2466        pts.sort()
     2467        return pts
    23542468
    23552469    ## def newform_eval(self, z, prec):
    23562470##         """
    class EllipticCurve_rational_field(Ellip 
    33843498        we return v = -1.
    33853499
    33863500        INPUT:
    3387             D (int) -- (deault: 0) Heegner discriminant; if 0, use the
     3501            D (int) -- (default: 0) Heegner discriminant; if 0, use the
    33883502                       first discriminant < -4 that satisfies the Heegner
    33893503                       hypothesis
    33903504            verbose (bool) -- (default: True)
    class EllipticCurve_rational_field(Ellip 
    35643678        Eq = ell_tate_curve.TateCurve(self,p)
    35653679        self._tate_curve[p] = Eq
    35663680        return Eq
     3681
     3682    def height(self, precision=53):
     3683        """
     3684        Returns the real height of this elliptic curve.
     3685        This is used in integral_points()
     3686
     3687
     3688        INPUT:
     3689
     3690            precision -- (integer: default 53) bit precision of the
     3691               field of real numbers in which the result will lie
     3692       
     3693        EXAMPLES:
     3694            sage: E=EllipticCurve('5077a1')
     3695            sage: E.height()               
     3696            17.4513334798896
     3697            sage: E.height(100)
     3698            17.451333479889612702508579399
     3699            sage: E=EllipticCurve([0,0,0,0,1])
     3700            sage: E.height()                 
     3701            1.38629436111989
     3702            sage: E=EllipticCurve([0,0,0,1,0])
     3703            sage: E.height()                 
     3704            7.45471994936400
     3705
     3706        """
     3707        R = RealField(precision)
     3708        c4 = self.c4()
     3709        c6 = self.c6()
     3710        j = self.j_invariant()
     3711        log_g2 = R((c4/12)).abs().log()
     3712        log_g3 = R((c6/216)).abs().log()
     3713
     3714        if j == 0:
     3715            h_j = R(1)
     3716        else:
     3717            h_j = max(log(R(abs(j.numerator()))), log(R(abs(j.denominator()))))
     3718        if (self.c4() != 0) and (self.c6() != 0):
     3719            h_gs = max(1, log_g2, log_g3)
     3720        elif c4 == 0:
     3721            if c6 == 0:
     3722                return max(1,h_j)
     3723            h_gs = max(1, log_g3)
     3724        else:
     3725            h_gs = max(1, log_g2)
     3726        return max(R(1),h_j, h_gs)
     3727
     3728    def antilogarithm(self, z, precision=100):
     3729        """
     3730        Returns the rational point (if any) associated to this complex
     3731        number; the inverse of the elliptic logarithm function.
     3732
     3733        INPUT:
     3734        z -- a complex number representing an element of
     3735             CC/L where L is the period lattice of the elliptic
     3736             curve
     3737
     3738        precision -- an integer specifying the precision (in bits) which
     3739                will be used for the computation
     3740
     3741        OUTPUT: The rational point which is the image of z under
     3742        the Weierstrass parametrization, if it exists and can
     3743        be determined from z with default precision.
     3744
     3745        NOTE: This uses the function ellztopoint from the pari library
     3746
     3747        TODO: Extend the wrapping of ellztopoint() to allow passing of
     3748        the precision parameter.
     3749           
     3750        """
     3751       
     3752        E_pari = self.pari_curve(prec)
     3753        try:
     3754            coords = E_pari.ellztopoint(z)
     3755            if len(coords) == 1:
     3756                return self(0)
     3757            return self([CC(xy).real().simplest_rational() for xy in coords])
     3758        except TypeError:
     3759            raise ValueError, "No rational point computable from z"
     3760
     3761    def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
     3762        """
     3763        Computes all integral points (up to sign) on this elliptic curve.
     3764       
     3765        INPUT:
     3766            mw_base -- list of EllipticCurvePoint generating the
     3767                       Mordell-Weil group of E
     3768                    (default: 'auto' - calls self.gens())
     3769            both_signs -- True/False (default False): if True the
     3770                       output contains both P and -P, otherwise only
     3771                       one of each pair.         
     3772            verbose -- True/False (default False): if True, some
     3773                       details of the computation are output
     3774                                   
     3775        OUTPUT:
     3776            A sorted list of all the integral points on E
     3777            (up to sign unless both_signs is True)
     3778           
     3779        NOTES: The complexity increases exponentially in the rank of
     3780            curve E.  The computation time (but not the output!)
     3781            depends on the Mordell-Weil basis.  If mw_base is given
     3782            but it not a basis for the Mordell-Weil group (modulo
     3783            torsion), integral points which are not in the subgroup
     3784            generated by the given points will almost certainly not be
     3785            listed.
     3786           
     3787        EXAMPLES:
     3788        A curve of rank 3 with no torsion points
     3789       
     3790            sage: E=EllipticCurve([0,0,1,-7,6])
     3791            sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
     3792            sage: a=E.integral_points([P1,P2,P3]); a 
     3793            [(-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)]
     3794           
     3795            sage: a = E.integral_points([P1,P2,P3], verbose=True)
     3796            Using mw_basis  [(2 : 0 : 1), (4 : 6 : 1), (114/49 : -720/343 : 1)]
     3797            e1,e2,e3:  -3.01243037259331 1.06582054769620 1.94660982489710
     3798            Minimal eigenvalue of height pairing matrix:  0.472730555831538
     3799            x-coords of points on compact component with  -3 <=x<= 1
     3800            set([0, -1, -3, -2, 1])
     3801            x-coords of points on non-compact component with  2 <=x<= 6
     3802            set([2, 3, 4])
     3803            starting search of remaining points using coefficient bound  6
     3804            x-coords of extra integral points:
     3805            set([2, 3, 4, 37, 406, 8, 11, 14, 816, 52, 21, 342, 93])
     3806            Total number of integral points: 18
     3807
     3808        It is also possible to not specify mw_base, but then the
     3809        Mordell-Weil basis must be computed;  this may take much longer
     3810       
     3811            sage: E=EllipticCurve([0,0,1,-7,6])
     3812            sage: a=E.integral_points(both_signs=True); a
     3813            [(-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)]
     3814           
     3815
     3816        An example with negative discriminant:
     3817            sage: EllipticCurve('900d1').integral_points()
     3818            [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)]
     3819       
     3820        Another example with rank 5 and no torsion points
     3821           
     3822            sage: E=EllipticCurve([-879984,319138704])
     3823            sage: P1=E.point((540,1188)); P2=E.point((576,1836))
     3824            sage: P3=E.point((468,3132)); P4=E.point((612,3132))
     3825            sage: P5=E.point((432,4428))
     3826            sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (400s!)
     3827            54
     3828           
     3829           
     3830        NOTES:
     3831            - This function uses the algorithm given in [Co1]
     3832            REFERENCES:
     3833                -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations
     3834                        GTM 239, Springer 2007
     3835        AUTHORS:
     3836            - Michael Mardaus (2008-07)
     3837            - Tobias Nagel (2008-07)
     3838            - John Cremona (2008-07)
     3839        """
     3840        #####################################################################
     3841        # INPUT CHECK #######################################################
     3842        if not self.is_integral():
     3843            raise ValueError, "integral_points() can only be called on an integral model"
     3844       
     3845        if mw_base=='auto':
     3846            mw_base = self.gens()
     3847            r = len(mw_base)
     3848        else:
     3849            try:
     3850                r = len(mw_base)
     3851            except TypeError:
     3852                raise TypeError, 'mw_base must be a list'
     3853            if not all([P.curve() is self for P in mw_base]):
     3854                raise ValueError, "points are not on the correct curve"
     3855               
     3856        tors_points = self.torsion_points()
     3857
     3858        # END INPUT-CHECK####################################################
     3859        #####################################################################
     3860               
     3861        #####################################################################
     3862        # INTERNAL FUNCTIONS ################################################
     3863   
     3864        ############################## begin ################################
     3865        def point_preprocessing(list):
     3866            #Transforms mw_base so that at most one point is on the
     3867            #compact component of the curve
     3868            Q = []
     3869            mem = -1
     3870            for i in range(0,len(list)):
     3871                if not list[i].is_on_identity_component(): # i.e. is on "egg"
     3872                    if mem == -1:
     3873                        mem = i
     3874                    else:
     3875                        Q.append(list[i]+list[mem])
     3876                        mem = i
     3877                else:
     3878                    Q.append(list[i])
     3879            if mem != -1: #add last point, which is not in egg, to Q
     3880                Q.append(2*list[mem])
     3881            return Q
     3882        ##############################  end  ################################
     3883        ############################## begin ################################
     3884        def modified_height(i):#computes modified height if base point i
     3885            return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2)
     3886        ##############################  end  ################################
     3887        ############################## begin ################################
     3888        def integral_x_coords_in_interval(xmin,xmax):
     3889            """
     3890            Returns the set of integers x with xmin<=x<=xmax which are
     3891            x-coordinates of points on this curve.
     3892            """
     3893            return set([x for x  in range(xmin,xmax) if self.is_x_coord(x)])
     3894        ##############################  end  ################################
     3895        ############################## begin ################################
     3896        def integral_points_with_bounded_mw_ceoffs():
     3897            r"""
     3898            Returns the set of integers x which are x-coordinates of
     3899            points on the curve which are linear combinations of the
     3900            generators (basis and torsion points) with coefficients
     3901            bounded by $H_q$.  The bound $H_q$ will be computed at
     3902            runtime.
     3903            """
     3904            from sage.misc.all import cartesian_product_iterator
     3905            from sage.groups.generic import multiples
     3906            xs=set()
     3907            N=H_q
     3908           
     3909            def use(P):
     3910                """
     3911                Helper function to record x-coord of a point if integral.
     3912                """
     3913                if not P.is_zero():
     3914                    xP = P[0]
     3915                    if xP.is_integral():
     3916                        xs.add(xP)
     3917                       
     3918            def use_t(R):
     3919                """
     3920                Helper function to record x-coords of a point +torsion if integral.
     3921                """
     3922                for T in tors_points:
     3923                    use(R+T)
     3924                       
     3925            # We use a naive method when the number of possibilities is small:
     3926
     3927            if r==1 or N<=10:
     3928                for P in multiples(mw_base[0],N+1):
     3929                    use_t(P)
     3930                return xs
     3931           
     3932            # Otherwise is is very very much faster to first compute
     3933            # the linear combinations over RR, and only compte them as
     3934            # rational points if they are approximately integral
     3935
     3936            def is_approx_integral(P):
     3937                return (abs(P[0]-P[0].round()))<0.001 and (abs(P[1]-P[1].round()))<0.001
     3938           
     3939            RR = RealField() #(100)
     3940            ER = self.change_ring(RR)
     3941            Rgens = [ER.lift_x(P[0]) for P in mw_base]
     3942            for i in range(r):
     3943                if abs(Rgens[i][1]-mw_base[i][1])>abs((-Rgens[i])[1]-mw_base[i][1]):
     3944                    Rgens[i] = -Rgens[i]
     3945            for ni in cartesian_product_iterator([range(-N,N+1) for i in range(r-1)]+[range(N+1)]):
     3946                RP=sum([ni[i]*Rgens[i] for i in range(r)],ER(0))
     3947                for T in tors_points:
     3948                    if is_approx_integral(RP+ER(T)):
     3949                        P = sum([ni[i]*mw_base[i] for i in range(r)],T)
     3950                        use(P)
     3951            return xs
     3952           
     3953            # Below this point we keep earlier experimental code which
     3954            # is slower when either N or r is large
     3955           
     3956#           if r==2:
     3957#               for P in multiples(mw_base[0],N+1):
     3958#                   for Q in multiples(mw_base[1],N+1):
     3959#                       for R in set([Q+P,Q-P]):
     3960#                           use_t(R)
     3961#               return xs
     3962#           
     3963#           if r==3:
     3964#               for P in multiples(mw_base[0],N+1):
     3965#                   for Q in multiples(mw_base[1],N+1):
     3966#                       PQ = [P+Q,P-Q]
     3967#                       PQ = set(PQ + [-R for R in PQ]) # {\pm P\pm Q}
     3968#                       for R in multiples(mw_base[2],N+1):
     3969#                           for S in PQ:
     3970#                               use_t(R+S)
     3971#               return xs
     3972#           
     3973#           # general rank
     3974#           mults=[list(multiples(mw_base[i],N+1)) for i in range(r)]
     3975#           for i in range(r-1):
     3976#               mults[i] = [-P for P in mults[i] if not P.is_zero()] + mults[i]
     3977#           for Pi in cartesian_product_iterator(mults):
     3978#               use_t(sum(Pi,self(0)))
     3979#           return xs
     3980#
     3981#  older, even slower  code:
     3982#
     3983#           for i in range(ceil(((2*H_q+1)**r)/2)):
     3984#               koeffs = Z(i).digits(base=2*H_q+1)
     3985#               koeffs = [0]*(r-len(koeffs)) + koeffs # pad with 0s
     3986#               P = sum([(koeffs[j]-H_q)*mw_base[j] for j in range(r)],self(0))
     3987#               for Q in tors_points: # P + torsion points (includes 0)
     3988#                   tmp = P + Q
     3989#                   if not tmp.is_zero():
     3990#                       x = tmp[0]
     3991#                       if x.is_integral():
     3992#                           xs.add(x)
     3993#           return xs
     3994        ##############################  end  #################################
     3995   
     3996        # END Internal functions #############################################
     3997        ######################################################################
     3998   
     3999        if (r == 0):
     4000            int_points = [P for P in tors_points if not P.is_zero()]
     4001            int_points = [P for P in int_points if P[0].is_integral()]
     4002            if not both_signs:
     4003                xlist = set([P[0] for P in int_points])
     4004                int_points = [self.lift_x(x) for x in xlist]
     4005            int_points.sort()
     4006            if verbose:
     4007                print 'Total number of integral points:',len(int_points)
     4008            return int_points
     4009
     4010
     4011        g2 = self.c4()/12
     4012        g3 = self.c6()/216
     4013        disc = self.discriminant()
     4014        j = self.j_invariant()
     4015        b2 = self.b2()
     4016       
     4017        Qx = rings.PolynomialRing(RationalField(),'x')
     4018        pol = Qx([-self.c6()/216,-self.c4()/12,0,4])
     4019        if disc > 0: # two real component -> 3 roots in RR
     4020            #on curve 897e4, only one root is found with default precision!
     4021            RR = R
     4022            prec = RR.precision()
     4023            ei = pol.roots(RR,multiplicities=False)
     4024            while len(ei)<3:
     4025                prec*=2
     4026                RR=RealField(prec)
     4027                ei = pol.roots(RR,multiplicities=False)
     4028            e1,e2,e3 = ei
     4029            if r >= 2: #preprocessing of mw_base only necessary if rank > 1
     4030                mw_base = point_preprocessing(mw_base) #at most one point in
     4031                                                       #E^{egg}, saved in P_egg
     4032
     4033        elif disc < 0: # one real component => 1 root in RR (=: e3),
     4034                       # 2 roots in C (e1,e2)
     4035            roots = pol.roots(C,multiplicities=False)
     4036            e3 = pol.roots(R,multiplicities=False)[0]
     4037            roots.remove(e3)
     4038            e1,e2 = roots
     4039 
     4040        e = R(1).exp()
     4041        pi = R(constants.pi)
     4042   
     4043        if verbose:
     4044            print "Using mw_basis ",mw_base
     4045            print "e1,e2,e3: ",e1,e2,e3
     4046           
     4047        # Algorithm presented in [Co1]
     4048        h_E = self.height()
     4049        w1, w2 = self.period_lattice().basis()
     4050        mu = R(disc).abs().log() / 6
     4051        if j!=0:
     4052            mu += max(R(1),R(j).abs().log()) / 6
     4053        if b2!=0:
     4054            mu += max(R(1),R(b2).abs().log())
     4055            mu += log(R(2))
     4056        else:
     4057            mu += 1
     4058       
     4059        c1 = (mu + 2.14).exp()
     4060        M = self.height_pairing_matrix(mw_base)
     4061        c2 = min(M.charpoly ().roots(multiplicities=False))
     4062        if verbose:
     4063            print "Minimal eigenvalue of height pairing matrix: ", c2
     4064
     4065        c3 = (w1**2)*R(b2).abs()/48 + 8
     4066        c5 = (c1*c3).sqrt()
     4067        c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
     4068       
     4069        mw_base_log = [] #contains \Phi(Q_i)
     4070        mod_h_list = []  #contains h_m(Q_i)
     4071        c9_help_list = []
     4072        for i in range(0,r):
     4073            mw_base_log.append(mw_base[i].elliptic_logarithm().abs())
     4074            mod_h_list.append(modified_height(i))
     4075            c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i])
     4076   
     4077        c8 = max(e*h_E,max(mod_h_list))
     4078        c9 = e/c7.sqrt() * min(c9_help_list)
     4079        n=r+1
     4080        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))
     4081   
     4082        top = 128 #arbitrary first upper bound
     4083        bottom = 0
     4084        log_c9=log(c9); log_c5=log(c5)
     4085        log_r_top = log(R(r*(10**top)))
     4086#        if verbose:
     4087#            print "[bottom,top] = ",[bottom,top]
     4088           
     4089        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):
     4090            #initial bound 'top' too small, upshift of search interval
     4091            bottom = top
     4092            top = 2*top
     4093        while top >= bottom: #binary-search like search for fitting exponent (bound)
     4094#            if verbose:
     4095#                print "[bottom,top] = ",[bottom,top]
     4096            bound = floor(bottom + (top - bottom)/2)
     4097            log_r_bound = log(R(r*(10**bound)))
     4098            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):
     4099                bottom = bound + 1
     4100            else:
     4101                top = bound - 1
     4102   
     4103        H_q = R(10)**bound
     4104        break_cond = 0 #at least one reduction step
     4105        #reduction via LLL
     4106        while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10%
     4107            c = R((H_q**n)*10)  #c has to be greater than H_q^n
     4108            M = matrix.MatrixSpace(Z,n)
     4109            m = M.identity_matrix()
     4110            for i in range(r):
     4111                m[i, r] = R(c*mw_base_log[i]).round()
     4112            m[r,r] = R(c*w1).round()
     4113   
     4114            #LLL - implemented in sage - operates on rows not on columns
     4115            m_LLL = m.LLL()
     4116            m_gram = m_LLL.gram_schmidt()[0]
     4117            b1_norm = R(m_LLL.row(0).norm())
     4118   
     4119            #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0
     4120            c1_LLL = -1
     4121            for i in range(n):
     4122                tmp = R(b1_norm/(m_gram.row(i).norm()))
     4123                if tmp > c1_LLL:
     4124                    c1_LLL = tmp
     4125   
     4126            if c1_LLL < 0:
     4127                raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base'
     4128           
     4129            d_L_0 = R(b1_norm**2 / c1_LLL)
     4130   
     4131            #Reducing of upper bound
     4132            Q = r * H_q**2
     4133            T = (1 + (3/2*r*H_q))/2
     4134            if d_L_0 < R(T**2+Q):
     4135                d_L_0 = 10*(T**2*Q)
     4136            low_bound = R((sqrt(d_L_0 - Q) - T)/c)
     4137   
     4138            #new bound according to low_bound and upper bound
     4139            #[c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
     4140            if low_bound != 0:
     4141                H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2)))
     4142                H_q_new = ceil(H_q_new)
     4143                if H_q_new == 1:
     4144                    break_cond = 1 # stops reduction
     4145                else:
     4146                    break_cond = R(H_q_new/H_q)
     4147                H_q = H_q_new
     4148            else:
     4149                break_cond = 1 # stops reduction, so using last H_q > 0
     4150            #END LLL-Reduction loop
     4151   
     4152        b2_12 = b2/12
     4153        if disc > 0:
     4154            ##Points in egg have X(P) between e1 and e2 [X(P)=x(P)+b2/12]:
     4155            x_int_points = integral_x_coords_in_interval(ceil(e1-b2_12), floor(e2-b2_12)+1)
     4156            if verbose:
     4157                print 'x-coords of points on compact component with ',ceil(e1-b2_12),'<=x<=',floor(e2-b2_12)
     4158                print x_int_points
     4159        else:
     4160            x_int_points = set()
     4161           
     4162        ##Points in noncompact component with X(P)< 2*max(|e1|,|e2|,|e3|) , espec. X(P)>=e3
     4163        x0 = ceil(e3-b2_12)
     4164        x1 = ceil(2*max(abs(e1),abs(e2),abs(e3)) - b2_12)
     4165        x_int_points2 = integral_x_coords_in_interval(x0, x1)
     4166        x_int_points = x_int_points.union(x_int_points2)
     4167        if verbose:
     4168            print 'x-coords of points on non-compact component with ',x0,'<=x<=',x1-1
     4169            print x_int_points2
     4170       
     4171        if verbose:
     4172            print 'starting search of remaining points using coefficient bound ',H_q
     4173        x_int_points3 = integral_points_with_bounded_mw_ceoffs()
     4174        x_int_points = x_int_points.union(x_int_points3)
     4175        if verbose:
     4176            print 'x-coords of extra integral points:'
     4177            print x_int_points3
     4178
     4179        # Now we have all the x-coordinates of integral points, and we
     4180        # construct the points, depending on the parameter both_signs:
     4181        if both_signs:
     4182            int_points = sum([self.lift_x(x,all=True) for x in x_int_points],[])
     4183        else:
     4184            int_points = [self.lift_x(x) for x in x_int_points]
     4185        int_points.sort()
     4186        if verbose:
     4187            print 'Total number of integral points:',len(int_points)
     4188        return int_points
    35674189
    35684190
    35694191def cremona_curves(conductors):
  • sage/schemes/elliptic_curves/rational_torsion.py

    diff -r 354bcfd4fd38 sage/schemes/elliptic_curves/rational_torsion.py
    a b class EllipticCurveTorsionSubgroup(group 
    4848
    4949    AUTHOR: Nick Alexander
    5050    """
    51     def __init__(self, E, flag=0):
     51    def __init__(self, E, algorithm="pari"):
    5252        self.__E = E
    5353
     54        flag = 0
     55        if algorithm=="doud":
     56            flag = 1
     57        if algorithm=="lutz_nagell":
     58            flag = 2
    5459        G = None
    5560        loop = 0
    5661        while G is None and loop < 3: