Ticket #1115: sage-trac1115.patch

File sage-trac1115.patch, 20.3 KB (added by cremona, 13 years ago)
  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1220268320 -3600
    # Node ID 34eba3d2cf2d23d5ead18c854bcc6148f4690bf9
    # Parent  d3af9775143a072f8bd0fd1babd7bb2777efa1d0
    #1115: enhanced precision handling for heights, regulators, analytic sha
    
    diff -r d3af9775143a -r 34eba3d2cf2d sage/schemes/elliptic_curves/ell_point.py
    a b  
    5555#*****************************************************************************
    5656
    5757from sage.structure.element import AdditiveGroupElement, RingElement
    58 
     58from sage.interfaces import gp
    5959import sage.plot.all as plot
    6060
    6161from sage.rings.padics.factory import Qp
     
    646646        ans.sort()
    647647        return ans
    648648
    649     def height(self):
     649    def height(self, precision=None):
    650650        """
    651651        The Neron-Tate canonical height of the point.
    652652
    653653        INPUT:
    654654            self -- a point on a curve over QQ
     655
     656            precision -- (int or None): the desired precision in bits;
     657                                        default real precision of None.
    655658
    656659        OUTPUT:
    657660            the rational number 0 or a nonzero real field element
     
    689692            Elliptic Curve defined by y^2 + 17*x*y - 120*y = x^3 - 60*x^2 over Rational Field
    690693            sage: E([30, -90]).height()
    691694            0
     695
     696            sage: E = EllipticCurve('389a1'); E
     697            Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field
     698            sage: [P,Q] = [E(-1,1),E(0,-1)]
     699            sage: P.height(precision=100)
     700            0.68666708330558658572355210295
     701            sage: (3*Q).height(precision=100)/Q.height(precision=100)
     702            9.0000000000000000000000000000
     703            sage: _.parent()
     704            Real Field with 100 bits of precision
    692705        """
    693706        if self.is_finite_order():
    694707            return rings.QQ(0)
    695         h = self.curve().pari_curve().ellheight([self[0], self[1]])
    696         return rings.RR(h)
     708
     709        if precision is None:
     710            precision = rings.RealField().precision()
     711        prec = (rings.RealField()(0.301029995663981)*precision).ceil() # decimal precision
     712
     713        # Note: we construct the pari_curve with the correct
     714        # precision; then it is not necessary to pass the precision on
     715        # to the height function.
     716        h = self.curve().pari_curve(prec).ellheight([self[0], self[1]])
     717        return rings.RealField(precision)(h)
    697718   
    698719    def is_on_identity_component(self):
    699720        """
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r d3af9775143a -r 34eba3d2cf2d sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    12601260    def gens(self, verbose=False, rank1_search=10,
    12611261             algorithm='mwrank_shell',
    12621262             only_use_mwrank=True,
    1263              proof = None):
     1263             proof = None,
     1264             use_database = True):
    12641265        """
    12651266        Compute and return generators for the Mordell-Weil group E(Q)
    12661267        *modulo* torsion.
     
    12901291                       first use more naive, natively implemented methods.
    12911292            proof -- bool or None (default None, see proof.elliptic_curve or
    12921293                       sage.structure.proof). 
    1293         OUTPUT:
    1294             generators -- List of generators for the Mordell-Weil group.
     1294            use_database -- bool (default True) if True, attempts to
     1295                       find curve and gens in the (optional) database
     1296        OUTPUT:
     1297            generators -- List of generators for the Mordell-Weil
     1298                          group modulo torsion.
    12951299
    12961300        IMPLEMENTATION: Uses Cremona's mwrank C library.
    12971301
     
    13021306
    13031307        A non-integral example:
    13041308            sage: E = EllipticCurve([-3/8,-2/3])
    1305             sage: E.gens()
     1309            sage: E.gens() # random (up to sign)
    13061310            [(10/9 : 29/54 : 1)]
    13071311
    13081312        A non-minimal example:
    13091313            sage: E = EllipticCurve('389a1')
    13101314            sage: E1 = E.change_weierstrass_model([1/20,0,0,0]); E1
    13111315            Elliptic Curve defined by y^2 + 8000*y = x^3 + 400*x^2 - 320000*x over Rational Field
    1312             sage: E1.gens()
     1316            sage: E1.gens() # random (if database not used)
    13131317            [(-400 : 8000 : 1), (0 : -8000 : 1)]
    13141318        """
    13151319        if proof is None:
     
    13171321            proof = get_flag(proof, "elliptic_curve")
    13181322        else:
    13191323            proof = bool(proof)
     1324
     1325        # If the gens are already cached, return them:
    13201326        try:
    13211327            return list(self.__gens[proof])  # return copy so not changed
    13221328        except KeyError:
    13231329            if proof is False and self.__gens.has_key(True):
    13241330                return self.__gens[True]
     1331
     1332        # If the optional extended database is installed and an
     1333        # isomorphic curve is in the database then its gens will be
     1334        # known:
     1335
     1336        if use_database:
     1337            try:
     1338                E = self.database_curve()
     1339                iso = E.isomorphism_to(self)
     1340                self.__gens[True] = [iso(P) for P in E.gens()]
     1341                return self.__gens[True]
     1342            except (RuntimeError, KeyError):  # curve or gens not in database
     1343                pass
     1344       
    13251345        if self.conductor() > 10**7:
    13261346            only_use_mwrank = True
    13271347
     
    13331353                if r == 0:
    13341354                    misc.verbose("Rank = 0, so done.")
    13351355                    self.__gens[True] = []
    1336                     self.__regulator[True] = R(1)
    13371356                    return self.__gens[True]
    13381357                if r == 1 and rank1_search:
    13391358                    misc.verbose("Rank = 1, so using direct search.")
     
    13441363                        G = [P for P in G if P.order() == oo]
    13451364                        if len(G) > 0:
    13461365                            misc.verbose("Direct search succeeded.")
    1347                             G, _, reg = self.saturation(G, verbose=verbose)
     1366                            G, _, _ = self.saturation(G, verbose=verbose)
    13481367                            misc.verbose("Computed saturation.")
    13491368                            self.__gens[True] = G
    1350                             self.__regulator[True] = reg
    13511369                            return self.__gens[True]
    13521370                        h += 2
    13531371                    misc.verbose("Direct search FAILED.")
     
    13621380            G = C.gens()
    13631381            if proof is True and C.certain() is False:
    13641382                del self.__mwrank_curve
    1365                 raise RuntimeError, "Unable to compute the rank, hence generators, with certainty (lower bound=%s).  This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + \
     1383                raise RuntimeError, "Unable to compute the rank, hence generators, with certainty (lower bound=%s, generators found=%s).  This could be because Sha(E/Q)[2] is nontrivial."%(C.rank(),G) + \
    13661384                      "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again."
    13671385            else:
    13681386                proof = C.certain()
     
    13951413                G.append(eval(X[k:j].replace(':',',')))
    13961414                X = X[j:]
    13971415                i = X.find('Generator ')
    1398             i = X.find('Regulator = ')
    1399             j = i + X[i:].find('\n')
    1400             self.__regulator[proof] = R(X[i+len('Regulator = '):j])
    14011416        ####
    14021417        self.__gens[proof] = [self.point(x, check=True) for x in G]
    14031418        self.__gens[proof].sort()
     
    14441459        """
    14451460        return len(self.gens(proof = proof))
    14461461
    1447     def regulator(self, use_database=True, verbose=None, proof=None):
     1462    def regulator(self, use_database=True, proof=None, precision=None):
    14481463        """
    14491464        Returns the regulator of this curve, which must be defined
    14501465        over Q.
    14511466
    14521467        INPUT:
    14531468            use_database -- bool (default: False), if True, try to
    1454                   look up the regulator in the Cremona database.
    1455             verbose -- (default: None), if specified changes the
    1456                   verbosity of mwrank computations.
     1469                  look up the generators in the Cremona database.
    14571470            proof -- bool or None (default: None, see proof.[tab] or
    14581471                       sage.structure.proof).  Note that results from
    14591472                       databases are considered proof = True
     1473            precision -- int or Noene (default: None): the precision
     1474                         in bits of the result (default real precision
     1475                         if None)
    14601476
    14611477        EXAMPLES:
    14621478            sage: E = EllipticCurve([0, 0, 1, -1, 0])
     
    14751491            sage: EllipticCurve([0, 0, 1, -79, 342]).regulator(proof=False)  # long time (seconds)
    14761492            14.790527570131...
    14771493        """
     1494        if precision is None:
     1495            RR = rings.RealField()
     1496            precision = RR.precision()
     1497        else:
     1498            RR = rings.RealField(precision)
     1499
    14781500        if proof is None:
    14791501            from sage.structure.proof.proof import get_flag
    14801502            proof = get_flag(proof, "elliptic_curve")
    14811503        else:
    14821504            proof = bool(proof)
    1483         try:
    1484             return self.__regulator[proof]
     1505
     1506        # We return a cached value if it exists and has sufficient precision:
     1507        try:
     1508            reg = self.__regulator[proof]
     1509            if reg.parent().precision() >= precision:
     1510                return RR(reg)
     1511            else: # Found regulator value but precision is too low
     1512                pass
    14851513        except KeyError:
    14861514            if proof is False and self.__regulator.has_key(True):
    1487                 return self.__regulator[True]
    1488         if use_database:
    1489             try:
    1490                 self.__regulator[True] = R(self.database_curve().db_extra[3])
    1491                 return self.__regulator[True]
    1492             except (AttributeError, RuntimeError):
    1493                 pass
    1494         G = self.gens(proof=proof)
    1495         try:  # in some cases self.gens() efficiently computes regulator.
    1496             return self.__regulator[proof]
    1497         except KeyError:
    1498             if proof is False and self.__regulator.has_key(True):
    1499                 return self.__regulator[True]
    1500         C = self.mwrank_curve()
    1501         reg = R(C.regulator())
    1502         if proof is True and not C.certain():
    1503             raise RuntimeError, "Unable to compute the rank, hence regulator, with certainty (lower bound=%s)."%C.rank()
    1504         proof = C.certain()
    1505         self.__regulator[proof] = reg
     1515                reg = self.__regulator[True]
     1516                if reg.parent().precision() >= precision:
     1517                    return RR(reg)
     1518                else: # Found regulator value but precision is too low
     1519                    pass
     1520       
     1521        # Next we find the gens, taking them from the database if they
     1522        # are there and use_database is True, else computing them:   
     1523
     1524        G = self.gens(proof=proof, use_database=use_database)
     1525
     1526        # Finally compute the regulator of the generators found:
     1527
     1528        self.__regulator[proof] = self.regulator_of_points(G,precision=precision)
    15061529        return self.__regulator[proof]
    15071530
    15081531    def height_pairing_matrix(self, points=None, precision=None):
     
    15161539                  curve, or (default) None, in which case self.gens()
    15171540                  will be used.
    15181541        precision -- number of bit of precision of result
    1519                  (default: default RealField)
     1542                 (default: None, for default RealField precision)
    15201543                 
    15211544        TODO: implement variable precision for heights of points
    15221545
     
    15291552            sage: EllipticCurve('11a').height_pairing_matrix()
    15301553            []
    15311554            sage: E=EllipticCurve('5077a1')
    1532             sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]])             
    1533             [  1.36857250535393  -1.30957670708658 -0.634867157837156]
    1534             [ -1.30957670708658   2.71735939281229   1.09981843056673]
    1535             [-0.634867157837156   1.09981843056673  0.668205165651928]
    1536 
     1555            sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]], precision=100)
     1556            [  1.3685725053539301576677189587  -1.3095767070865762526921116660 -0.63486715783715585992297292250]
     1557            [ -1.3095767070865762526921116660   2.7173593928122929952451158897   1.0998184305667293436670206574]
     1558            [-0.63486715783715585992297292250   1.0998184305667293436670206574  0.66820516565192789038007958879]
    15371559        """
    15381560        if points is None:
    15391561            points = self.gens()
     
    15561578                mat[k,j]=mat[j,k]
    15571579        return mat
    15581580
     1581    def regulator_of_points(self, points=[], precision=None):
     1582        """
     1583        Returns the regulator of the given points on this curve.
     1584
     1585        INPUT:
     1586
     1587        points -- a list of points on this curve (default: empty list)
     1588        precision -- number of bit of precision of result
     1589                 (default: None, for default RealField precision)
     1590                 
     1591        TODO: implement variable precision for heights of points
     1592
     1593        EXAMPLES:
     1594            sage: E = EllipticCurve('37a1')                 
     1595            sage: P = E(0,0)
     1596            sage: Q = E(1,0)
     1597            sage: E.regulator_of_points([P,Q])
     1598            0.000000000000000
     1599            sage: 2*P==Q
     1600            True
     1601           
     1602            sage: E = EllipticCurve('5077a1')
     1603            sage: points = [E.lift_x(x) for x in [-2,-7/4,1]]
     1604            sage: E.regulator_of_points(points)             
     1605            0.417143558758384
     1606            sage: E.regulator_of_points(points,precision=100)
     1607            0.41714355875838396981711954462
     1608
     1609            sage: E = EllipticCurve('389a')                   
     1610            sage: E.regulator_of_points()             
     1611            1.00000000000000
     1612            sage: points = [P,Q] = [E(-1,1),E(0,-1)]         
     1613            sage: E.regulator_of_points(points)               
     1614            0.152460177943144
     1615            sage: E.regulator_of_points(points, precision=100)
     1616            0.15246017794314375162432475705
     1617            sage: E.regulator_of_points(points, precision=200)
     1618            0.15246017794314375162432475704945582324372707748663081784028
     1619            sage: E.regulator_of_points(points, precision=300)
     1620            0.152460177943143751624324757049455823243727077486630817840280980046053225683562463604114816
     1621        """
     1622        for P in points:
     1623            assert P.curve() == self
     1624
     1625        if precision is None:
     1626            RR = rings.RealField()
     1627            precision = RR.precision()
     1628        else:
     1629            RR = rings.RealField(precision)
     1630        r = len(points)
     1631        M = matrix.MatrixSpace(RR, r)
     1632        mat = M()
     1633        for j in range(r):
     1634            mat[j,j] = points[j].height(precision=precision)
     1635        for j in range(r):
     1636            for k in range(j+1,r):             
     1637                mat[j,k]=((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2
     1638                mat[k,j]=mat[j,k]
     1639        return mat.det()
     1640
    15591641
    15601642    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False):
    15611643        """
     
    15751657            saturation (list) -- points that form a basis for the saturation
    15761658            index (int) -- the index of the group generated by points
    15771659                           in their saturation
    1578             regulator (float) -- regulator of saturated points.
     1660            regulator (real with default precision) -- regulator of saturated points.
    15791661
    15801662        IMPLEMENTATION:
    15811663            Uses Cremona's mwrank package.  With max_prime=0, we call
     
    19752057        try:
    19762058            return self.__tamagawa_product
    19772059        except AttributeError:
    1978             self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python()
     2060            self.__tamagawa_product = Integer(self.pari_mincurve().ellglobalred()[2].python())
    19792061            return self.__tamagawa_product
    19802062
    19812063    def real_components(self):
     
    36853767        Returns the real height of this elliptic curve.
    36863768        This is used in integral_points()
    36873769
    3688 
    3689         INPUT:
    3690 
     3770        INPUT:
    36913771            precision -- (integer: default 53) bit precision of the
    36923772               field of real numbers in which the result will lie
    36933773       
  • sage/schemes/elliptic_curves/sha.py

    diff -r d3af9775143a -r 34eba3d2cf2d sage/schemes/elliptic_curves/sha.py
    a b  
    11from sage.structure.sage_object import SageObject
    22from sage.rings.all import (
    33    Integer,
     4    RealField,
    45    RationalField,
    56    RIF)
    67from sage.misc.functional import log
     
    2021    ########################################################################
    2122
    2223    def an_numerical(self, prec = 53,
    23                          use_database=False, proof=None):
     24                         use_database=True, proof=None):
    2425        """
    2526        Return the numerical analytic order of Sha, which is
    2627        a floating point number in all cases.
    2728
    2829        INPUT:
    29             prec -- integer (default: 53) bits precision -- just used
    30                     for the L-series computation; not for regulator, etc.
    31             use_database -- whether the rank and regulator should
     30            prec -- integer (default: 53) bits precision -- used
     31                    for the L-series computation, period, regulator, etc.
     32            use_database -- whether the rank and generators should
    3233                    be looked up in the database if possible.
    3334            proof -- bool or None (default: None, see proof.[tab] or
    3435                           sage.structure.proof) proof option passed
     
    3637
    3738        NOTE: See also the an() command, which will return a
    3839        provably correct integer when the rank is 0 or 1.
     40
     41        WARNING: If the curve's generators are not known, computing
     42        them may be very time-consuming.  Also, computation of the
     43        L-series derivative will be time-consuming for large rank and
     44        large conductor, and the computation time for this may
     45        increase substantially at greater precision.  However, use of
     46        very low precision less than about 10 can cause the underlying
     47        pari library functions to fail.
    3948
    4049        EXAMPLES:
    4150            sage: EllipticCurve('11a').sha().an_numerical()
     
    5362            sage: EllipticCurve([1, -1, 0, -79, 289]).sha().an_numerical()   # long time
    5463            1.00000000000000
    5564
    56     A rank 5 curve:
    57             sage: EllipticCurve([0, 0, 1, -79, 342]).sha().an_numerical(prec=4, proof=False)          # long time -- about 30 seconds.
    58             1.0               
     65        A rank 5 curve:
     66            sage: EllipticCurve([0, 0, 1, -79, 342]).sha().an_numerical(prec=10, proof=False)          # long time -- about 30 seconds.
     67            1.0
     68
     69            # See trac #1115
     70            sage: sha=EllipticCurve('37a1').sha()                       
     71            sage: [sha.an_numerical(prec) for prec in xrange(30,100,10)]
     72            [1.0000000,
     73            1.0000000000,
     74            1.0000000000000,
     75            1.0000000000000000,
     76            1.0000000000000000000,
     77            1.0000000000000000000000,
     78            1.0000000000000000000000000]
    5979        """
     80        RR = RealField(prec)
     81        prec2 = prec+2
     82        RR2 = RealField(prec2)
    6083        try:
    61             return self.__an_numerical
     84            an = self.__an_numerical
     85            if an.parent().precision() >= prec:
     86                return RR(an)
     87            else: # cached precision too low
     88                pass
    6289        except AttributeError:
    6390            pass
    6491        r = Integer(self.E.rank(use_database=use_database, proof=proof))
    65         L = self.E.lseries().dokchitser(prec=prec)
    66         Lr= L.derivative(1,r)
    67         Om = self.E.period_lattice().omega()
    68         Reg = self.E.regulator(use_database=use_database, proof=proof)
     92        L = self.E.lseries().dokchitser(prec=prec2)
     93        Lr= RR2(L.derivative(1,r))  # L.derivative() returns a Complex
     94        Om = RR2(self.E.period_lattice().omega(prec2))
     95        Reg = self.E.regulator(use_database=use_database, proof=proof, precision=prec2)
    6996        T = self.E.torsion_order()
    7097        cp = self.E.tamagawa_product()
    71         Sha = Lr*T*T/(r.factorial()*Om*cp*Reg)
     98        Sha = RR((Lr*T*T)/(r.factorial()*Om*cp*Reg))
    7299        self.__an_numerical = Sha
    73100        return Sha
    74101