Ticket #1115: 1115sha_prec.patch
File 1115sha_prec.patch, 20.7 KB (added by , 13 years ago) 


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 56 56 #***************************************************************************** 57 57 58 58 from sage.structure.element import AdditiveGroupElement, RingElement 59 59 from sage.interfaces import gp 60 60 import sage.plot.all as plot 61 61 62 62 from sage.rings.padics.factory import Qp … … 976 976 gxdd = gxd.derivative() 977 977 return ( e(gxd(self[0])) > 0 and e(gxdd(self[0])) > 0) 978 978 979 def height(self ):979 def height(self, precision=None): 980 980 """ 981 981 The NeronTate canonical height of the point. 982 982 … … 985 985 986 986 INPUT: 987 987 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) 988 992 989 993 OUTPUT: 990 994 the rational number 0 or a nonzero real field element … … 1022 1026 Elliptic Curve defined by y^2 + 17*x*y  120*y = x^3  60*x^2 over Rational Field 1023 1027 sage: E([30, 90]).height() 1024 1028 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 """ 1026 1040 if self.has_finite_order(): 1027 1041 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. 1028 1050 try: 1029 h = self.curve().pari_curve( ).ellheight([self[0], self[1]])1030 return rings.R R(h)1051 h = self.curve().pari_curve(prec).ellheight([self[0], self[1]]) 1052 return rings.RealField(precision)(h) 1031 1053 except AttributeError: "canonical height not yet implemented over general number fields." 1032 1054 1033 1055 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 1260 1260 def gens(self, verbose=False, rank1_search=10, 1261 1261 algorithm='mwrank_shell', 1262 1262 only_use_mwrank=True, 1263 proof = None): 1263 proof = None, 1264 use_database = True): 1264 1265 """ 1265 1266 Compute and return generators for the MordellWeil group E(Q) 1266 1267 *modulo* torsion. … … 1290 1291 first use more naive, natively implemented methods. 1291 1292 proof  bool or None (default None, see proof.elliptic_curve or 1292 1293 sage.structure.proof). 1293 OUTPUT: 1294 generators  List of generators for the MordellWeil 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 MordellWeil 1298 group modulo torsion. 1295 1299 1296 1300 IMPLEMENTATION: Uses Cremona's mwrank C library. 1297 1301 … … 1302 1306 1303 1307 A nonintegral example: 1304 1308 sage: E = EllipticCurve([3/8,2/3]) 1305 sage: E.gens() 1309 sage: E.gens() # random (up to sign) 1306 1310 [(10/9 : 29/54 : 1)] 1307 1311 1308 1312 A nonminimal example: 1309 1313 sage: E = EllipticCurve('389a1') 1310 1314 sage: E1 = E.change_weierstrass_model([1/20,0,0,0]); E1 1311 1315 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) 1313 1317 [(400 : 8000 : 1), (0 : 8000 : 1)] 1314 1318 """ 1315 1319 if proof is None: … … 1317 1321 proof = get_flag(proof, "elliptic_curve") 1318 1322 else: 1319 1323 proof = bool(proof) 1324 1325 # If the gens are already cached, return them: 1320 1326 try: 1321 1327 return list(self.__gens[proof]) # return copy so not changed 1322 1328 except KeyError: 1323 1329 if proof is False and self.__gens.has_key(True): 1324 1330 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 1325 1345 if self.conductor() > 10**7: 1326 1346 only_use_mwrank = True 1327 1347 … … 1333 1353 if r == 0: 1334 1354 misc.verbose("Rank = 0, so done.") 1335 1355 self.__gens[True] = [] 1336 self.__regulator[True] = R(1)1337 1356 return self.__gens[True] 1338 1357 if r == 1 and rank1_search: 1339 1358 misc.verbose("Rank = 1, so using direct search.") … … 1344 1363 G = [P for P in G if P.order() == oo] 1345 1364 if len(G) > 0: 1346 1365 misc.verbose("Direct search succeeded.") 1347 G, _, reg= self.saturation(G, verbose=verbose)1366 G, _, _ = self.saturation(G, verbose=verbose) 1348 1367 misc.verbose("Computed saturation.") 1349 1368 self.__gens[True] = G 1350 self.__regulator[True] = reg1351 1369 return self.__gens[True] 1352 1370 h += 2 1353 1371 misc.verbose("Direct search FAILED.") … … 1362 1380 G = C.gens() 1363 1381 if proof is True and C.certain() is False: 1364 1382 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) + \ 1366 1384 "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again." 1367 1385 else: 1368 1386 proof = C.certain() … … 1395 1413 G.append(eval(X[k:j].replace(':',','))) 1396 1414 X = X[j:] 1397 1415 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])1401 1416 #### 1402 1417 self.__gens[proof] = [self.point(x, check=True) for x in G] 1403 1418 self.__gens[proof].sort() … … 1444 1459 """ 1445 1460 return len(self.gens(proof = proof)) 1446 1461 1447 def regulator(self, use_database=True, verbose=None, proof=None):1462 def regulator(self, use_database=True, proof=None, precision=None): 1448 1463 """ 1449 1464 Returns the regulator of this curve, which must be defined 1450 1465 over Q. 1451 1466 1452 1467 INPUT: 1453 1468 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. 1457 1470 proof  bool or None (default: None, see proof.[tab] or 1458 1471 sage.structure.proof). Note that results from 1459 1472 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) 1460 1476 1461 1477 EXAMPLES: 1462 1478 sage: E = EllipticCurve([0, 0, 1, 1, 0]) … … 1475 1491 sage: EllipticCurve([0, 0, 1, 79, 342]).regulator(proof=False) # long time (seconds) 1476 1492 14.790527570131... 1477 1493 """ 1494 if precision is None: 1495 RR = rings.RealField() 1496 precision = RR.precision() 1497 else: 1498 RR = rings.RealField(precision) 1499 1478 1500 if proof is None: 1479 1501 from sage.structure.proof.proof import get_flag 1480 1502 proof = get_flag(proof, "elliptic_curve") 1481 1503 else: 1482 1504 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 1485 1513 except KeyError: 1486 1514 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) 1506 1529 return self.__regulator[proof] 1507 1530 1508 1531 def height_pairing_matrix(self, points=None, precision=None): … … 1515 1538 points  either a list of points, which must be on this 1516 1539 curve, or (default) None, in which case self.gens() 1517 1540 will be used. 1518 precision  number of bit of precision of result1519 (default: default RealField)1541 precision  number of bits of precision of result 1542 (default: None, for default RealField precision) 1520 1543 1521 1544 TODO: implement variable precision for heights of points 1522 1545 … … 1529 1552 sage: EllipticCurve('11a').height_pairing_matrix() 1530 1553 [] 1531 1554 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] 1537 1559 """ 1538 1560 if points is None: 1539 1561 points = self.gens() … … 1556 1578 mat[k,j]=mat[j,k] 1557 1579 return mat 1558 1580 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 1559 1642 1560 1643 def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False): 1561 1644 """ … … 1575 1658 saturation (list)  points that form a basis for the saturation 1576 1659 index (int)  the index of the group generated by points 1577 1660 in their saturation 1578 regulator ( float)  regulator of saturated points.1661 regulator (real with default precision)  regulator of saturated points. 1579 1662 1580 1663 IMPLEMENTATION: 1581 1664 Uses Cremona's mwrank package. With max_prime=0, we call … … 1975 2058 try: 1976 2059 return self.__tamagawa_product 1977 2060 except AttributeError: 1978 self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python()2061 self.__tamagawa_product = Integer(self.pari_mincurve().ellglobalred()[2].python()) 1979 2062 return self.__tamagawa_product 1980 2063 1981 2064 def real_components(self): … … 3782 3865 Returns the real height of this elliptic curve. 3783 3866 This is used in integral_points() 3784 3867 3785 3786 INPUT: 3787 3868 INPUT: 3788 3869 precision  (integer: default 53) bit precision of the 3789 3870 field of real numbers in which the result will lie 3790 3871 
sage/schemes/elliptic_curves/sha.py
diff r 6e4ac02aec85 r ab31aa3ef3c6 sage/schemes/elliptic_curves/sha.py
a b 1 1 from sage.structure.sage_object import SageObject 2 2 from sage.rings.all import ( 3 3 Integer, 4 RealField, 4 5 RationalField, 5 6 RIF) 6 7 from sage.misc.functional import log … … 20 21 ######################################################################## 21 22 22 23 def an_numerical(self, prec = 53, 23 use_database= False, proof=None):24 use_database=True, proof=None): 24 25 """ 25 26 Return the numerical analytic order of Sha, which is 26 27 a floating point number in all cases. 27 28 28 29 INPUT: 29 prec  integer (default: 53) bits precision  justused30 for the Lseries computation ; not forregulator, etc.31 use_database  whether the rank and regulatorshould30 prec  integer (default: 53) bits precision  used 31 for the Lseries computation, period, regulator, etc. 32 use_database  whether the rank and generators should 32 33 be looked up in the database if possible. 33 34 proof  bool or None (default: None, see proof.[tab] or 34 35 sage.structure.proof) proof option passed … … 36 37 37 38 NOTE: See also the an() command, which will return a 38 39 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 timeconsuming. Also, computation of the 43 Lseries derivative will be timeconsuming 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. 39 48 40 49 EXAMPLES: 41 50 sage: EllipticCurve('11a').sha().an_numerical() … … 53 62 sage: EllipticCurve([1, 1, 0, 79, 289]).sha().an_numerical() # long time 54 63 1.00000000000000 55 64 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] 59 79 """ 80 RR = RealField(prec) 81 prec2 = prec+2 82 RR2 = RealField(prec2) 60 83 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 62 89 except AttributeError: 63 90 pass 64 91 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) 69 96 T = self.E.torsion_order() 70 97 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)) 72 99 self.__an_numerical = Sha 73 100 return Sha 74 101