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 #***************************************************************************** from sage.structure.element import AdditiveGroupElement, RingElement from sage.interfaces import gp import sage.plot.all as plot from sage.rings.padics.factory import Qp gxdd = gxd.derivative() return ( e(gxd(self[0])) > 0 and e(gxdd(self[0])) > 0) def height(self): def height(self, precision=None): """ The Neron-Tate canonical height of the point. INPUT: self -- a point on a curve over QQ precision -- int or None (default: None): the precision in bits of the result (default real precision if None) OUTPUT: the rational number 0 or a nonzero real field element Elliptic Curve defined by y^2 + 17*x*y - 120*y = x^3 - 60*x^2 over Rational Field sage: E([30, -90]).height() 0 """ sage: E = EllipticCurve('389a1'); E Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field sage: [P,Q] = [E(-1,1),E(0,-1)] sage: P.height(precision=100) 0.68666708330558658572355210295 sage: (3*Q).height(precision=100)/Q.height(precision=100) 9.0000000000000000000000000000 sage: _.parent() Real Field with 100 bits of precision """ if self.has_finite_order(): return rings.QQ(0) if precision is None: precision = rings.RealField().precision() prec = (rings.RealField()(0.301029995663981)*precision).ceil() # decimal precision # Note: we construct the pari_curve with the correct # precision; then it is not necessary to pass the precision on # to the height function. try: h = self.curve().pari_curve().ellheight([self[0], self[1]]) return rings.RR(h) h = self.curve().pari_curve(prec).ellheight([self[0], self[1]]) return rings.RealField(precision)(h) except AttributeError: "canonical height not yet implemented over general number fields." 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 def gens(self, verbose=False, rank1_search=10, algorithm='mwrank_shell', only_use_mwrank=True, proof = None): proof = None, use_database = True): """ Compute and return generators for the Mordell-Weil group E(Q) *modulo* torsion. first use more naive, natively implemented methods. proof -- bool or None (default None, see proof.elliptic_curve or sage.structure.proof). OUTPUT: generators -- List of generators for the Mordell-Weil group. use_database -- bool (default True) if True, attempts to find curve and gens in the (optional) database OUTPUT: generators -- List of generators for the Mordell-Weil group modulo torsion. IMPLEMENTATION: Uses Cremona's mwrank C library. A non-integral example: sage: E = EllipticCurve([-3/8,-2/3]) sage: E.gens() sage: E.gens() # random (up to sign) [(10/9 : 29/54 : 1)] A non-minimal example: sage: E = EllipticCurve('389a1') sage: E1 = E.change_weierstrass_model([1/20,0,0,0]); E1 Elliptic Curve defined by y^2 + 8000*y = x^3 + 400*x^2 - 320000*x over Rational Field sage: E1.gens() sage: E1.gens() # random (if database not used) [(-400 : 8000 : 1), (0 : -8000 : 1)] """ if proof is None: proof = get_flag(proof, "elliptic_curve") else: proof = bool(proof) # If the gens are already cached, return them: try: return list(self.__gens[proof])  # return copy so not changed except KeyError: if proof is False and self.__gens.has_key(True): return self.__gens[True] # If the optional extended database is installed and an # isomorphic curve is in the database then its gens will be # known: if use_database: try: E = self.database_curve() iso = E.isomorphism_to(self) self.__gens[True] = [iso(P) for P in E.gens(use_database=False)] return self.__gens[True] except (RuntimeError, KeyError):  # curve or gens not in database pass if self.conductor() > 10**7: only_use_mwrank = True if r == 0: misc.verbose("Rank = 0, so done.") self.__gens[True] = [] self.__regulator[True] = R(1) return self.__gens[True] if r == 1 and rank1_search: misc.verbose("Rank = 1, so using direct search.") G = [P for P in G if P.order() == oo] if len(G) > 0: misc.verbose("Direct search succeeded.") G, _, reg = self.saturation(G, verbose=verbose) G, _, _ = self.saturation(G, verbose=verbose) misc.verbose("Computed saturation.") self.__gens[True] = G self.__regulator[True] = reg return self.__gens[True] h += 2 misc.verbose("Direct search FAILED.") G = C.gens() if proof is True and C.certain() is False: del self.__mwrank_curve 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() + \ 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) + \ "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again." else: proof = C.certain() G.append(eval(X[k:j].replace(':',','))) X = X[j:] i = X.find('Generator ') i = X.find('Regulator = ') j = i + X[i:].find('\n') self.__regulator[proof] = R(X[i+len('Regulator = '):j]) #### self.__gens[proof] = [self.point(x, check=True) for x in G] self.__gens[proof].sort() """ return len(self.gens(proof = proof)) def regulator(self, use_database=True, verbose=None, proof=None): def regulator(self, use_database=True, proof=None, precision=None): """ Returns the regulator of this curve, which must be defined over Q. INPUT: use_database -- bool (default: False), if True, try to look up the regulator in the Cremona database. verbose -- (default: None), if specified changes the verbosity of mwrank computations. look up the generators in the Cremona database. proof -- bool or None (default: None, see proof.[tab] or sage.structure.proof).  Note that results from databases are considered proof = True precision -- int or None (default: None): the precision in bits of the result (default real precision if None) EXAMPLES: sage: E = EllipticCurve([0, 0, 1, -1, 0]) sage: EllipticCurve([0, 0, 1, -79, 342]).regulator(proof=False)  # long time (seconds) 14.790527570131... """ if precision is None: RR = rings.RealField() precision = RR.precision() else: RR = rings.RealField(precision) if proof is None: from sage.structure.proof.proof import get_flag proof = get_flag(proof, "elliptic_curve") else: proof = bool(proof) try: return self.__regulator[proof] # We return a cached value if it exists and has sufficient precision: try: reg = self.__regulator[proof] if reg.parent().precision() >= precision: return RR(reg) else: # Found regulator value but precision is too low pass except KeyError: if proof is False and self.__regulator.has_key(True): return self.__regulator[True] if use_database: try: self.__regulator[True] = R(self.database_curve().db_extra[3]) return self.__regulator[True] except (AttributeError, RuntimeError): pass G = self.gens(proof=proof) try:  # in some cases self.gens() efficiently computes regulator. return self.__regulator[proof] except KeyError: if proof is False and self.__regulator.has_key(True): return self.__regulator[True] C = self.mwrank_curve() reg = R(C.regulator()) if proof is True and not C.certain(): raise RuntimeError, "Unable to compute the rank, hence regulator, with certainty (lower bound=%s)."%C.rank() proof = C.certain() self.__regulator[proof] = reg reg = self.__regulator[True] if reg.parent().precision() >= precision: return RR(reg) else: # Found regulator value but precision is too low pass # Next we find the gens, taking them from the database if they # are there and use_database is True, else computing them: G = self.gens(proof=proof, use_database=use_database) # Finally compute the regulator of the generators found: self.__regulator[proof] = self.regulator_of_points(G,precision=precision) return self.__regulator[proof] def height_pairing_matrix(self, points=None, precision=None): points -- either a list of points, which must be on this curve, or (default) None, in which case self.gens() will be used. precision -- number of bit of precision of result (default: default RealField) precision -- number of bits of precision of result (default: None, for default RealField precision) TODO: implement variable precision for heights of points sage: EllipticCurve('11a').height_pairing_matrix() [] sage: E=EllipticCurve('5077a1') sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]]) [  1.36857250535393  -1.30957670708658 -0.634867157837156] [ -1.30957670708658   2.71735939281229   1.09981843056673] [-0.634867157837156   1.09981843056673  0.668205165651928] sage: E.height_pairing_matrix([E.lift_x(x) for x in [-2,-7/4,1]], precision=100) [  1.3685725053539301576677189587  -1.3095767070865762526921116660 -0.63486715783715585992297292250] [ -1.3095767070865762526921116660   2.7173593928122929952451158897   1.0998184305667293436670206574] [-0.63486715783715585992297292250   1.0998184305667293436670206574  0.66820516565192789038007958879] """ if points is None: points = self.gens() mat[k,j]=mat[j,k] return mat def regulator_of_points(self, points=[], precision=None): """ Returns the regulator of the given points on this curve. INPUT: points -- a list of points on this curve (default: empty list) precision -- int or None (default: None): the precision in bits of the result (default real precision if None) TODO: implement variable precision for heights of points EXAMPLES: sage: E = EllipticCurve('37a1') sage: P = E(0,0) sage: Q = E(1,0) sage: E.regulator_of_points([P,Q]) 0.000000000000000 sage: 2*P==Q True sage: E = EllipticCurve('5077a1') sage: points = [E.lift_x(x) for x in [-2,-7/4,1]] sage: E.regulator_of_points(points) 0.417143558758384 sage: E.regulator_of_points(points,precision=100) 0.41714355875838396981711954462 sage: E = EllipticCurve('389a') sage: E.regulator_of_points() 1.00000000000000 sage: points = [P,Q] = [E(-1,1),E(0,-1)] sage: E.regulator_of_points(points) 0.152460177943144 sage: E.regulator_of_points(points, precision=100) 0.15246017794314375162432475705 sage: E.regulator_of_points(points, precision=200) 0.15246017794314375162432475704945582324372707748663081784028 sage: E.regulator_of_points(points, precision=300) 0.152460177943143751624324757049455823243727077486630817840280980046053225683562463604114816 """ for P in points: assert P.curve() == self if precision is None: RR = rings.RealField() precision = RR.precision() else: RR = rings.RealField(precision) r = len(points) M = matrix.MatrixSpace(RR, r) mat = M() for j in range(r): mat[j,j] = points[j].height(precision=precision) for j in range(r): for k in range(j+1,r): mat[j,k]=((points[j]+points[k]).height(precision=precision) - mat[j,j] - mat[k,k])/2 mat[k,j]=mat[j,k] return mat.det() def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False): """ saturation (list) -- points that form a basis for the saturation index (int) -- the index of the group generated by points in their saturation regulator (float) -- regulator of saturated points. regulator (real with default precision) -- regulator of saturated points. IMPLEMENTATION: Uses Cremona's mwrank package.  With max_prime=0, we call try: return self.__tamagawa_product except AttributeError: self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python() self.__tamagawa_product = Integer(self.pari_mincurve().ellglobalred()[2].python()) return self.__tamagawa_product def real_components(self): Returns the real height of this elliptic curve. This is used in integral_points() INPUT: INPUT: precision -- (integer: default 53) bit precision of the field of real numbers in which the result will lie
• sage/schemes/elliptic_curves/sha.py

`diff -r 6e4ac02aec85 -r ab31aa3ef3c6 sage/schemes/elliptic_curves/sha.py`
 a from sage.structure.sage_object import SageObject from sage.rings.all import ( Integer, RealField, RationalField, RIF) from sage.misc.functional import log ######################################################################## def an_numerical(self, prec = 53, use_database=False, proof=None): use_database=True, proof=None): """ Return the numerical analytic order of Sha, which is a floating point number in all cases. INPUT: prec -- integer (default: 53) bits precision -- just used for the L-series computation; not for regulator, etc. use_database -- whether the rank and regulator should prec -- integer (default: 53) bits precision -- used for the L-series computation, period, regulator, etc. use_database -- whether the rank and generators should be looked up in the database if possible. proof -- bool or None (default: None, see proof.[tab] or sage.structure.proof) proof option passed NOTE: See also the an() command, which will return a provably correct integer when the rank is 0 or 1. WARNING: If the curve's generators are not known, computing them may be very time-consuming.  Also, computation of the L-series derivative will be time-consuming for large rank and large conductor, and the computation time for this may increase substantially at greater precision.  However, use of very low precision less than about 10 can cause the underlying pari library functions to fail. EXAMPLES: sage: EllipticCurve('11a').sha().an_numerical() sage: EllipticCurve([1, -1, 0, -79, 289]).sha().an_numerical()   # long time 1.00000000000000 A rank 5 curve: sage: EllipticCurve([0, 0, 1, -79, 342]).sha().an_numerical(prec=4, proof=False)          # long time -- about 30 seconds. 1.0 A rank 5 curve: sage: EllipticCurve([0, 0, 1, -79, 342]).sha().an_numerical(prec=10, proof=False)          # long time -- about 30 seconds. 1.0 # See trac #1115 sage: sha=EllipticCurve('37a1').sha() sage: [sha.an_numerical(prec) for prec in xrange(30,100,10)] [1.0000000, 1.0000000000, 1.0000000000000, 1.0000000000000000, 1.0000000000000000000, 1.0000000000000000000000, 1.0000000000000000000000000] """ RR = RealField(prec) prec2 = prec+2 RR2 = RealField(prec2) try: return self.__an_numerical an = self.__an_numerical if an.parent().precision() >= prec: return RR(an) else: # cached precision too low pass except AttributeError: pass r = Integer(self.E.rank(use_database=use_database, proof=proof)) L = self.E.lseries().dokchitser(prec=prec) Lr= L.derivative(1,r) Om = self.E.period_lattice().omega() Reg = self.E.regulator(use_database=use_database, proof=proof) L = self.E.lseries().dokchitser(prec=prec2) Lr= RR2(L.derivative(1,r))  # L.derivative() returns a Complex Om = RR2(self.E.period_lattice().omega(prec2)) Reg = self.E.regulator(use_database=use_database, proof=proof, precision=prec2) T = self.E.torsion_order() cp = self.E.tamagawa_product() Sha = Lr*T*T/(r.factorial()*Om*cp*Reg) Sha = RR((Lr*T*T)/(r.factorial()*Om*cp*Reg)) self.__an_numerical = Sha return Sha