Ticket #3674: sage-trac3674c.2.patch

File sage-trac3674c.2.patch, 19.2 KB (added by cremona, 9 years ago)

replaces earlier sage-trac3674c.patch

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1217382417 14400
    # Node ID 4bb3dbd8745648726c2682b39070b1d78b0c7ddb
    # Parent  f9aaca6611fb37297e8e4d730e7cd1205f638ceb
    trac#3674:  more work on integral points & related functions
    
    diff -r f9aaca6611fb -r 4bb3dbd87456 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 
    867866                y = lam*(x-x3)-y-a1*x3-a3
    868867                x = x3
    869868            c = (x - e3).sqrt()
    870             while (a-b)/a > 0.5**precision:
     869            while (a - b)/a > 0.5**(precision-1):
    871870                a,b,c = (a + b)/2, (a*b).sqrt(), (c + (c**2 + b**2 - a**2).sqrt())/2
    872871
    873872            z = (a/c).arcsin()/a
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r f9aaca6611fb -r 4bb3dbd87456 sage/schemes/elliptic_curves/ell_rational_field.py
    a b AUTHORS: 
    1212   -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files.
    1313   -- John Cremona (2008-01)
    1414   -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points
     15   -- John Cremona (2008-07): further work on integral_points
    1516"""
    1617
    1718#*****************************************************************************
    import sage.databases.cremona 
    5657import sage.databases.cremona
    5758from   sage.libs.pari.all import pari
    5859import sage.functions.transcendental as transcendental
    59 import math
    60 from sage.calculus.calculus import sqrt, arcsin, floor, ceil
     60from sage.calculus.calculus import sqrt, floor, ceil
    6161import sage.libs.mwrank.all as mwrank
    6262import constructor
    6363from sage.interfaces.all import gp
    import ell_tate_curve 
    8686import ell_tate_curve
    8787
    8888factor = arith.factor
    89 exp = math.exp
    9089mul = misc.mul
    9190next_prime = arith.next_prime
    9291kronecker_symbol = arith.kronecker_symbol
    class EllipticCurve_rational_field(Ellip 
    23522351            self.__torsion_order = self.torsion_subgroup().order()
    23532352            return self.__torsion_order
    23542353
    2355     def torsion_subgroup(self, flag=0):
     2354    def torsion_subgroup(self, algorithm="pari"):
    23562355        """
    23572356        Returns the torsion subgroup of this elliptic curve.
    23582357
    23592358        INPUT:
    2360             flag -- (default: 0)  chooses PARI algorithm:
    2361               flag = 0: uses Doud algorithm
    2362               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
    23632363
    23642364        OUTPUT:
    23652365            The EllipticCurveTorsionSubgroup instance associated to this elliptic curve.
    class EllipticCurve_rational_field(Ellip 
    23902390        try:
    23912391            return self.__torsion_subgroup
    23922392        except AttributeError:
    2393             self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag)
     2393            self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, algorithm)
    23942394            self.__torsion_order = self.__torsion_subgroup.order()
    23952395            return self.__torsion_subgroup
    23962396
    2397     def torsion_points(self, flag=0):
     2397    def torsion_points(self, algorithm="pari"):
    23982398        """
    23992399        Returns the torsion points of this elliptic curve as a sorted list.
    24002400
    24012401        INPUT:
    2402             flag -- (default: 0)  chooses PARI algorithm:
    2403               flag = 0: uses Doud algorithm
    2404               flag = 1: uses Lutz-Nagell algorithm
     2402            algorithm -- string:
     2403                 "pari"  -- (default) use the pari library
     2404                 "doud" -- use Doud's algorithm
     2405                 "lutz_nagell" -- use the Lutz-Nagell theorem
    24052406
    24062407        OUTPUT:
    24072408            A list of all the torsion points on this elliptic curve.
    class EllipticCurve_rational_field(Ellip 
    24132414            [(0 : 1 : 0), (8 : -19 : 1), (8 : 18 : 1)]
    24142415
    24152416            sage: E=EllipticCurve([-1386747,368636886])   
    2416             sage: E.torsion_subgroup()
     2417            sage: T=E.torsion_subgroup(); T
    24172418            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
    24182423            sage: E.torsion_points() 
    24192424            [(-1293 : 0 : 1),
    24202425            (-933 : -29160 : 1),
    class EllipticCurve_rational_field(Ellip 
    24342439            (8787 : 816480 : 1)]
    24352440        """
    24362441        multiples = sage.groups.generic.multiples
    2437         gens = self.torsion_subgroup(flag).gens()
     2442        gens = self.torsion_subgroup(algorithm).gens()
    24382443        ngens = len(gens)
    24392444        if ngens == 0:
    24402445            return [self(0)]
    class EllipticCurve_rational_field(Ellip 
    34783483        we return v = -1.
    34793484
    34803485        INPUT:
    3481             D (int) -- (deault: 0) Heegner discriminant; if 0, use the
     3486            D (int) -- (default: 0) Heegner discriminant; if 0, use the
    34823487                       first discriminant < -4 that satisfies the Heegner
    34833488                       hypothesis
    34843489            verbose (bool) -- (default: True)
    class EllipticCurve_rational_field(Ellip 
    36603665        return Eq
    36613666
    36623667    def height(self, precision=53):
    3663         """Returns real height of this elliptic curve
     3668        """
     3669        Returns the real height of this elliptic curve.
    36643670        This is used in integral_points()
    36653671
    36663672
    3667         INPUT: precision (integer) bit precision of the field of real
    3668                numbers in which the result will lie (default: 53 as in
    3669                RealField())
     3673        INPUT:
     3674
     3675            precision -- (integer: default 53) bit precision of the
     3676               field of real numbers in which the result will lie
    36703677       
    36713678        EXAMPLES:
    36723679            sage: E=EllipticCurve('5077a1')
    class EllipticCurve_rational_field(Ellip 
    36823689            7.45471994936400
    36833690
    36843691        """
    3685         #Computes height of curve 'self' according to Co1
    36863692        R = RealField(precision)
    36873693        c4 = self.c4()
    36883694        c6 = self.c6()
    class EllipticCurve_rational_field(Ellip 
    37043710            h_gs = max(1, log_g2)
    37053711        return max(R(1),h_j, h_gs)
    37063712
    3707     def integral_points(self, mw_base='auto', both_signs=False):
    3708         """
    3709         Computes all integral points (up to sign) on the elliptic
    3710         curve E which has Mordell-Weil basis mw_base.
    3711        
    3712         INPUT:
    3713             self -- EllipticCurve_Rational_Field
     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:
    37143751            mw_base -- list of EllipticCurvePoint generating the
    37153752                       Mordell-Weil group of E
    37163753                    (default: 'auto' - calls self.gens())
    3717 
    37183754            both_signs -- True/False (default False): if True the
    37193755                       output contains both P and -P, otherwise only
    3720                        one of ecah pair.
    3721            
    3722         OUTPUT:
    3723             set of all integral points on E
     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
    37243762            (up to sign unless both_signs is True)
    37253763           
    3726         HINTS:
    3727             - The complexity increases exponentially in the rank of curve E.
    3728             - It can help if you try another Mordell-Weil base, because the
    3729             computation time depends on this, too.   
     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.
    37303771           
    37313772        EXAMPLES:
    37323773        A curve of rank 3 with no torsion points
    class EllipticCurve_rational_field(Ellip 
    37363777            sage: a=E.integral_points([P1,P2,P3]); a 
    37373778            [(-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)]
    37383779           
    3739         It is also possible to not specify mw_base and tors_base but
    3740         then the Mordell-Weil basis must be found by other functions
    3741         which will takes longer
     3780            sage: a = E.integral_points([P1,P2,P3], verbose=True)
     3781            Points on compact component:
     3782            [(-3 : -1 : 1), (-2 : -4 : 1), (-1 : -4 : 1), (0 : -3 : 1), (1 : -1 : 1)]
     3783            Points on non-compact component with x<= 7
     3784            [(2 : -1 : 1), (3 : -4 : 1), (4 : -7 : 1)]
     3785            starting search of remaining points using coefficient bound  9
     3786            Extra integral points:
     3787            [(342 : 6324 : 1), (37 : -225 : 1), (8 : 21 : 1), (52 : 374 : 1), (406 : -8181 : 1), (14 : -52 : 1), (816 : 23309 : 1), (93 : 896 : 1), (11 : -36 : 1), (3 : 3 : 1), (4 : -7 : 1), (21 : 95 : 1), (2 : -1 : 1)]
     3788            Total number of integral points: 18
     3789
     3790        It is also possible to not specify mw_base, but then the
     3791        Mordell-Weil basis must be computed;  this may take much longer
    37423792       
    37433793            sage: E=EllipticCurve([0,0,1,-7,6])
    37443794            sage: a=E.integral_points(both_signs=True); a
    class EllipticCurve_rational_field(Ellip 
    37713821        """
    37723822        #####################################################################
    37733823        # INPUT CHECK #######################################################
    3774         try:
    3775             assert self.is_integral()
    3776         except:
     3824        if not self.is_integral():
    37773825            raise ValueError, "integral_points() can only be called on an integral model"
    37783826       
    37793827        if mw_base=='auto':
    class EllipticCurve_rational_field(Ellip 
    37843832                r = len(mw_base)
    37853833            except TypeError:
    37863834                raise TypeError, 'mw_base must be a list'
    3787             for P in mw_base:
    3788                 try:
    3789                     assert P.curve() is self
    3790                 except:
    3791                     raise ValueError, "points are not on the correct curve"
     3835            if not all([P.curve() is self for P in mw_base]):
     3836                raise ValueError, "points are not on the correct curve"
    37923837               
    37933838        tors_base = self.torsion_subgroup().gens()
    37943839        len_tors = len(tors_base)
    class EllipticCurve_rational_field(Ellip 
    38343879        ##############################  end  ################################
    38353880        ############################## begin ################################
    38363881        def search_remaining_points():
    3837             """Returns list of integral points on curve E written as
    3838                linear combination of n times the mordell-weil base
    3839                points and torsion points (n is bounded by H_q, which
    3840                will be computed at runtime)
     3882            r"""
     3883            Returns a list of integral points on the curve which are
     3884            linear combinations of the generators (basis and torsion
     3885            points) with coefficients bounded by $H_q$.  The bound
     3886            $H_q$ will be computed at runtime.
    38413887            """
    3842             OP=self(0)
    38433888            pts=[]
    38443889            for i in range(1,((2*H_q+1)**r)/2):
    38453890                koeffs = Z(i).digits(base=2*H_q+1)
    class EllipticCurve_rational_field(Ellip 
    38493894                    P = P + (koeffs[index]-H_q)*mw_base[index]
    38503895                for Q in tors_points: # P + torsion points (includes 0)
    38513896                    P = P + Q
    3852                     if P!=OP and P[0].is_integral():
     3897                    if not P.is_zero() and P[0].is_integral():
    38533898                        pts.append(P)
    38543899                        if both_signs:
    38553900                            Q = -P
    class EllipticCurve_rational_field(Ellip 
    38683913            pts.sort()
    38693914            return pts
    38703915
    3871         e = exp(R(1))
     3916        e = R(1).exp()
    38723917        pi = R(constants.pi)
    38733918   
    38743919        g2 = self.c4()/12
    class EllipticCurve_rational_field(Ellip 
    38893934            e3 = pol.roots(R,multiplicities=False)[0]
    38903935            roots.remove(e3)
    38913936            e1,e2 = roots
    3892                
     3937
     3938        if verbose:
     3939            print "Using mw_basis ",mw_base
     3940            print "e1,e2,e3: ",e1,e2,e3
     3941           
    38933942        # Algorithm presented in [Co1]
    38943943        h_E = self.height()
    38953944        w1, w2 = self.period_lattice().basis()
    class EllipticCurve_rational_field(Ellip 
    39003949            mu += max(R(1),R(b2).abs().log())
    39013950            mu += log(R(2))
    39023951       
    3903         c1 = exp(mu + 2.14)
     3952        c1 = (mu + 2.14).exp()
    39043953        M = self.height_pairing_matrix(mw_base)
    39053954        c2 = min(M.charpoly ().roots(multiplicities=False))
     3955        if verbose:
     3956            print "Minimal eigenvalue of height pairing matrix: ", c2
     3957
    39063958        c3 = (w1**2)*R(b2).abs()/48 + 8
    39073959        c5 = (c1*c3).sqrt()
    39083960        c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
    class EllipticCurve_rational_field(Ellip 
    39243976        bottom = 0
    39253977        log_c9=log(c9); log_c5=log(c5)
    39263978        log_r_top = log(R(r*(10**top)))
    3927    
     3979        if verbose:
     3980            print "[bottom,top] = ",[bottom,top]
     3981           
    39283982        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):
    39293983            #initial bound 'top' too small, upshift of search interval
    39303984            bottom = top
    39313985            top = 2*top
    39323986        while top >= bottom: #binary-search like search for fitting exponent (bound)
     3987            if verbose:
     3988                print "[bottom,top] = ",[bottom,top]
    39333989            bound = floor(bottom + (top - bottom)/2)
    39343990            log_r_bound = log(R(r*(10**bound)))
    39353991            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):
    class EllipticCurve_rational_field(Ellip 
    39724028                d_L_0 = 10*(T**2*Q)
    39734029            low_bound = R((sqrt(d_L_0 - Q) - T)/c)
    39744030   
    3975             #new bound according to low_bound and upper bound [c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
     4031            #new bound according to low_bound and upper bound
     4032            #[c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
    39764033            H_q_new = R(sqrt(log(low_bound/c5)/(-c2/2)))
    39774034            H_q_new = ceil(H_q_new)
    39784035            break_cond = R(H_q_new/H_q)
    class EllipticCurve_rational_field(Ellip 
    39804037            #END LLL-Reduction loop
    39814038   
    39824039   
    3983         int_points = []
    3984         ##for a more detailed view how the function works uncomment
    3985         ##all print - statements below
    39864040        if disc > 0:
    3987             ##Points in egg have e2>=x>=e1
    3988             int_points += search_points(ceil(e1), floor(e2)+1)
    3989             #print 'Points on compact component \n',list(int_points)
     4041            ##Points in egg have x between e1 and e2:
     4042            int_points = search_points(ceil(e1), floor(e2)+1)
     4043            if verbose:
     4044                print 'Points on compact component:'
     4045                print int_points
     4046        else:
     4047            int_points = []
    39904048           
    39914049        ##Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3
    3992         int_points += search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3))))
    3993         #print 'Points on compact component and non-compact component part 1 \n', list(int_points)
    3994        
    3995         #print 'starting search of remainig points using bound ',H_q
    3996         int_points += search_remaining_points()
    3997         #print 'Integral points:\n',list(int_points)
    3998         #print 'Total amount of points:',len(int_points)
     4050        x0 = ceil(e3)
     4051        x1 = ceil(2*max(abs(e1),abs(e2),abs(e3)))
     4052        int_points2 = search_points(x0, x1)
     4053        int_points += int_points2
     4054        if verbose:
     4055            print 'Points on non-compact component with x<=',x1
     4056            print int_points2
     4057       
     4058        if verbose:
     4059            print 'starting search of remaining points using coefficient bound ',H_q
     4060        int_points3 = search_remaining_points()
     4061        int_points += int_points3
     4062        if verbose:
     4063            print 'Extra integral points:'
     4064            print int_points3
    39994065        if both_signs:
    40004066            int_points = list(set(int_points)) # remove duplicates
    40014067        else:
    40024068            xlist = set([P[0] for P in int_points])
    40034069            int_points = [self.lift_x(x) for x in xlist]
    40044070        int_points.sort()
     4071        if verbose:
     4072            print 'Total number of integral points:',len(int_points)
    40054073        return int_points
    40064074
    40074075
  • sage/schemes/elliptic_curves/rational_torsion.py

    diff -r f9aaca6611fb -r 4bb3dbd87456 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: