Ticket #3674: sage-trac3674c.2.patch

File sage-trac3674c.2.patch, 19.2 kB (added by cremona, 4 months ago)

replaces earlier sage-trac3674c.patch

  • a/sage/schemes/elliptic_curves/ell_point.py

    old new  
    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 
    …… 
    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 

  • a/sage/schemes/elliptic_curves/ell_rational_field.py

    old new  
    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#***************************************************************************** 
    …… 
    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 
    …… 
    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 
    …… 
    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. 
    …… 
    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. 
    …… 
    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), 
    …… 
    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)] 
    …… 
    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) 
    …… 
    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') 
    …… 
    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() 
    …… 
    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  
    …… 
    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 
    …… 
    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': 
    …… 
    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) 
    …… 
    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) 
    …… 
    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 
    …… 
    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 
    …… 
    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() 
    …… 
    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())) 
    …… 
    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): 
    …… 
    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) 
    …… 
    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 

  • a/sage/schemes/elliptic_curves/rational_torsion.py

    old new  
    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: