Ticket #1115: 1115-sha_prec.patch

File 1115-sha_prec.patch, 20.7 KB (added by AlexGhitza, 13 years ago)

apply only this patch, after #3377

  • sage/schemes/elliptic_curves/ell_point.py

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1220415322 -36000
    # Node ID ab31aa3ef3c6c3c9f879d58785fd79161e24f9d4
    # Parent  6e4ac02aec856889fb387a17d6a3c4392414f1e5
    trac 1115: elliptic_curves improvements to gens, regulator, sha
    
    diff -r 6e4ac02aec85 -r ab31aa3ef3c6 sage/schemes/elliptic_curves/ell_point.py
    a b  
    5656#*****************************************************************************
    5757
    5858from sage.structure.element import AdditiveGroupElement, RingElement
    59 
     59from sage.interfaces import gp
    6060import sage.plot.all as plot
    6161
    6262from sage.rings.padics.factory import Qp
     
    976976        gxdd = gxd.derivative()
    977977        return ( e(gxd(self[0])) > 0 and e(gxdd(self[0])) > 0)
    978978
    979     def height(self):
     979    def height(self, precision=None):
    980980        """
    981981        The Neron-Tate canonical height of the point.
    982982
     
    985985
    986986        INPUT:
    987987            self -- a point on a curve over QQ
     988
     989            precision -- int or None (default: None): the precision
     990                         in bits of the result (default real precision
     991                         if None)
    988992
    989993        OUTPUT:
    990994            the rational number 0 or a nonzero real field element
     
    10221026            Elliptic Curve defined by y^2 + 17*x*y - 120*y = x^3 - 60*x^2 over Rational Field
    10231027            sage: E([30, -90]).height()
    10241028            0
    1025         """
     1029
     1030            sage: E = EllipticCurve('389a1'); E
     1031            Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field
     1032            sage: [P,Q] = [E(-1,1),E(0,-1)]
     1033            sage: P.height(precision=100)
     1034            0.68666708330558658572355210295
     1035            sage: (3*Q).height(precision=100)/Q.height(precision=100)
     1036            9.0000000000000000000000000000
     1037            sage: _.parent()
     1038            Real Field with 100 bits of precision
     1039            """
    10261040        if self.has_finite_order():
    10271041            return rings.QQ(0)
     1042 
     1043        if precision is None:
     1044            precision = rings.RealField().precision()
     1045        prec = (rings.RealField()(0.301029995663981)*precision).ceil() # decimal precision
     1046
     1047        # Note: we construct the pari_curve with the correct
     1048        # precision; then it is not necessary to pass the precision on
     1049        # to the height function.
    10281050        try:
    1029             h = self.curve().pari_curve().ellheight([self[0], self[1]])
    1030             return rings.RR(h)
     1051            h = self.curve().pari_curve(prec).ellheight([self[0], self[1]])
     1052            return rings.RealField(precision)(h)
    10311053        except AttributeError: "canonical height not yet implemented over general number fields."
    10321054   
    10331055    def elliptic_logarithm(self, embedding=None, precision=100):
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r 6e4ac02aec85 -r ab31aa3ef3c6 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(use_database=False)]
     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 None (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):
     
    15151538        points -- either a list of points, which must be on this
    15161539                  curve, or (default) None, in which case self.gens()
    15171540                  will be used.
    1518         precision -- number of bit of precision of result
    1519                  (default: default RealField)
     1541        precision -- number of bits of precision of result
     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 -- int or None (default: None): the precision
     1589                         in bits of the result (default real precision
     1590                         if None)
     1591                 
     1592        TODO: implement variable precision for heights of points
     1593
     1594        EXAMPLES:
     1595            sage: E = EllipticCurve('37a1')                 
     1596            sage: P = E(0,0)
     1597            sage: Q = E(1,0)
     1598            sage: E.regulator_of_points([P,Q])
     1599            0.000000000000000
     1600            sage: 2*P==Q
     1601            True
     1602           
     1603            sage: E = EllipticCurve('5077a1')
     1604            sage: points = [E.lift_x(x) for x in [-2,-7/4,1]]
     1605            sage: E.regulator_of_points(points)             
     1606            0.417143558758384
     1607            sage: E.regulator_of_points(points,precision=100)
     1608            0.41714355875838396981711954462
     1609
     1610            sage: E = EllipticCurve('389a')                   
     1611            sage: E.regulator_of_points()             
     1612            1.00000000000000
     1613            sage: points = [P,Q] = [E(-1,1),E(0,-1)]         
     1614            sage: E.regulator_of_points(points)               
     1615            0.152460177943144
     1616            sage: E.regulator_of_points(points, precision=100)
     1617            0.15246017794314375162432475705
     1618            sage: E.regulator_of_points(points, precision=200)
     1619            0.15246017794314375162432475704945582324372707748663081784028
     1620            sage: E.regulator_of_points(points, precision=300)
     1621            0.152460177943143751624324757049455823243727077486630817840280980046053225683562463604114816
     1622        """
     1623        for P in points:
     1624            assert P.curve() == self
     1625
     1626        if precision is None:
     1627            RR = rings.RealField()
     1628            precision = RR.precision()
     1629        else:
     1630            RR = rings.RealField(precision)
     1631        r = len(points)
     1632        M = matrix.MatrixSpace(RR, r)
     1633        mat = M()
     1634        for j in range(r):
     1635            mat[j,j] = points[j].height(precision=precision)
     1636        for j in range(r):
     1637            for k in range(j+1,r):             
     1638                mat[j,k]=((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2
     1639                mat[k,j]=mat[j,k]
     1640        return mat.det()
     1641
    15591642
    15601643    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False):
    15611644        """
     
    15751658            saturation (list) -- points that form a basis for the saturation
    15761659            index (int) -- the index of the group generated by points
    15771660                           in their saturation
    1578             regulator (float) -- regulator of saturated points.
     1661            regulator (real with default precision) -- regulator of saturated points.
    15791662
    15801663        IMPLEMENTATION:
    15811664            Uses Cremona's mwrank package.  With max_prime=0, we call
     
    19752058        try:
    19762059            return self.__tamagawa_product
    19772060        except AttributeError:
    1978             self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python()
     2061            self.__tamagawa_product = Integer(self.pari_mincurve().ellglobalred()[2].python())
    19792062            return self.__tamagawa_product
    19802063
    19812064    def real_components(self):
     
    37823865        Returns the real height of this elliptic curve.
    37833866        This is used in integral_points()
    37843867
    3785 
    3786         INPUT:
    3787 
     3868        INPUT:
    37883869            precision -- (integer: default 53) bit precision of the
    37893870               field of real numbers in which the result will lie
    37903871       
  • sage/schemes/elliptic_curves/sha.py

    diff -r 6e4ac02aec85 -r ab31aa3ef3c6 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