Ticket #3883: sage-trac3883b.patch

File sage-trac3883b.patch, 39.9 kB (added by cremona, 4 months ago)

Apply after the previous patch

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

    old new  
    12431243 
    12441244############################################################# 
    12451245# 
    1246 # Explanation of the different division (also known as torsion) 
    1247 # polynomial functions in Sage 
     1246# Explanation of the division (also known as torsion) polynomial 
     1247# functions in Sage. 
    12481248# 
     1249# The main user function division_polynomial() (also aliased as 
     1250# torsion_polynomial()) is used to compute polynomials whose roots 
     1251# determine the $m$-torsion points on the curve.  Three options are 
     1252# available, which effect the result when $m$ is even and also the 
     1253# parent ring of the returned value.  The function can return either a 
     1254# polynomial or the evaluation of that polynomial at a point, 
     1255# depending on the input.  Values are cached. 
    12491256# 
    1250 # For odd $m$, the three functions pseudo_torsion_polynomial(m), 
    1251 # torsion_polynomial(m), full_division_polynomial(m) all return the 
    1252 # same polynomial of degree $(m^2-1)/2$ in $x$, though in the latter 
    1253 # case the returned polynomial is an element of the bivariate 
    1254 # polynomial ring in $x$ and $y$.  As a function on the curve, this 
    1255 # polynomial has $m^2-1$ zeros, namely the non-zero points $P$ such 
    1256 # that $m*P=0$.  (These degrees and multiplicities need adjusting in 
    1257 # finite characteristic $p$ when $p$ divides $m$.) 
    1258  
    1259 # For even $m$, the three functions are all different.  (1) 
    1260 # pseudo_torsion_polynomial(m) returns a polynomial in $x$ of degree 
    1261 # $(m^2-4)/2$ whose zeros (as a function on the curve) are the $m^2-4$ 
    1262 # points $P$ which satisfy $m*P=0$ but $2*P\not=0$.  (2) 
    1263 # torsion_polynomial(m) includes an extra factor of degree $3$ which 
    1264 # is zero at the $2$-torsion points, so as a function on the curve it 
    1265 # has a double zero at these points, while as a polynomial in $x$ it 
    1266 # has $(m^2+2)/2$ distinct roots which are the distinct 
    1267 # $x$-coordinates of the nonzero points $P$ satisfying $m*P=0$.  (3) 
    1268 # full_division_polynomial(m), which is a bivariate polynomial, has 
    1269 # instead an extra factor $2*y+a1*x+a3$, so as a function on the curve 
    1270 # it again has precisely $m^2-1$ zeros, namely the non-zero points $P$ 
    1271 # which satisfy $m*P=0$. 
    1272  
     1257# The options are controlled by the value of the parameter 
     1258# two_torsion_multiplicity, which may be 0, 1 or 2.  If it is 0 or 2, 
     1259# then a univariate polynomial will be returned (or evaluated at the 
     1260# parameter x if x is not None).  This is the polynomial whose roots 
     1261# are the values of $x(P)$ at the nonzero points $P$ where $m*P=0$ 
     1262# (when two_torsion_multiplicity==2), or the points where $m*P=0$ but 
     1263# $2*P\not=0$ (when two_torsion_multiplicity==0). 
     1264# 
     1265# If two_torsion_multiplicity==1, then a bivariate polynomial is 
     1266# returned, which (as a function on the curve) has a simple zero at 
     1267# each nonzero point $P$ such that $m*P=0$.  When $m$ is odd this is a 
     1268# polynomial in $x$ alone, but is still returned as an element of a 
     1269# polynomial ring in two variables; when $m$ is even it has a factor 
     1270# $2y+a_1x+a_3$.  In thise case if the parameter x is not None then it 
     1271# should be a tuple of length 2, or a point P on the curve, and the 
     1272# returned value is the value of the bivariate polynomial at this 
     1273# point. 
     1274# 
    12731275# Comparison with Magma: Magma's function DivisionPolynomial(E,m) 
    1274 # returns a triple of univariate polynomials.  For odd $m$ these are 
    1275 # f,f,1 where f is the same as both Sage's 
    1276 # E.pseudo_torsion_polynomial(m) and E.torsion_polynomial(m).  For 
    1277 # even $m$ they are f,g,h where f=E.torsion_polynomial(m), 
    1278 # g=E.pseudo_torsion_polynomial(m), and h=f/g.  Magma has no function 
    1279 # equivalent to Sage's E.full_division_polynomial(m) except for the 
    1280 # function TwoTorsionPolynomial(E) which is the same as 
    1281 # E.full_division_polynomial(2), i.e. 2*y+a1*x+a3. 
    1282  
    1283 # All these polynomials belong to the polynomial rings K[x] or K[x,y], 
    1284 # where K is the base field, and not to either the affine coordinate 
    1285 # ring of the curve or the function field. 
    1286  
    1287 # The pseudo_torsion_polynomial() function allows the user to supply a 
    1288 # value of 'x' which need not be the generator of a polynomial ring, 
    1289 # allowing for fast evaluation of the polynomials.  This feature is 
    1290 # used, for example, in the division_points() function in 
    1291 # ell_point.py.  There are implications for the caching of these 
    1292 # polynomials: see the docstrings of the individual functions for more 
    1293 # about this. 
     1276# returns a triple of univariate polynomials $f,g,h$ where $f$ is 
     1277# \code{E.division_polynomial(m,two_torsion_multiplicity=2)}, $g$ is 
     1278# \code{E.division_polynomial(m,two_torsion_multiplicity=0)} and $h$ 
     1279# is the quotient, so that $h=1$ when $m$ is odd. 
    12941280 
    12951281############################################################# 
    12961282 
    1297     def pseudo_torsion_polynomial(self, n, x=None, cache=None): 
     1283    def division_polynomial_0(self, n, x=None, cache=None): 
    12981284         r""" 
    1299          Returns the n-th torsion polynomial (division polynomial), without 
    1300          the 2-torsion factor if n is even, as a polynomial in $x$. 
     1285         Returns the $n$-th torsion (division) polynomial , without 
     1286         the 2-torsion factor if $n$ is even, as a polynomial in $x$. 
    13011287 
    13021288         These are the polynomials $g_n$ defined in Mazur/Tate (``The p-adic 
    13031289         sigma function''), but with the sign flipped for even $n$, so that 
    13041290         the leading coefficient is always positive. 
    13051291 
    1306          The full torsion polynomials may be recovered as follows 
    1307          (see the function full_division_polynomial()): 
    1308          \begin{itemize} 
    1309          \item $\psi_n = g_n$ for odd $n$. 
    1310          \item $\psi_n = (2y + a_1 x + a_3) g_n$ for even $n$. 
    1311          \end{itemize} 
    1312  
    1313          Note that the $g_n$'s are always polynomials in $x$, whereas the 
    1314          $\psi_n$'s require the appearance of a $y$. 
     1292         NOTE: This function is intended for internal use; users 
     1293         should use division_polynomial(). 
    13151294 
    13161295         SEE ALSO: 
    1317              -- torsion_polynomial() 
    13181296             -- multiple_x_numerator() 
    13191297             -- multiple_x_denominator() 
    1320              -- full_division_polynomial() 
     1298             -- division_polynomial() 
    13211299 
    13221300         INPUT: 
    13231301             n -- positive integer, or the special values -1 and -2 which 
    …… 
    13321310                  element of a ring. 
    13331311             cache -- optional dictionary, with integer keys. If the key m 
    13341312                  is in cache, then cache[m] is assumed to be the value of 
    1335                   pseudo_torsion_polynomial(m) for the supplied x. New entries 
     1313                  division_polynomial_0(m) for the supplied x. New entries 
    13361314                  will be added to the cache as they are computed. 
    13371315 
    13381316         ALGORITHM: 
    13391317             -- Recursion described in Mazur/Tate. The recursive formulae are 
    13401318             evaluated $O((log n)^2)$ times. 
    13411319 
    1342          TODO: 
    1343              -- for better unity of code, it might be good to make the regular 
    1344              torsion_polynomial() function use this as a subroutine.  At 
    1345              present they are independent. 
    1346  
    13471320         AUTHORS: 
    1348              -- David Harvey (2006-09-24) 
     1321             -- David Harvey (2006-09-24): initial version 
     1322             -- John Cremona (2008-08-26): unified division polynomial code 
    13491323 
    13501324         EXAMPLES: 
    13511325            sage: E = EllipticCurve("37a") 
    1352             sage: E.pseudo_torsion_polynomial(1) 
     1326            sage: E.division_polynomial_0(1) 
    13531327            1 
    1354             sage: E.pseudo_torsion_polynomial(2) 
     1328            sage: E.division_polynomial_0(2) 
    13551329            1 
    1356             sage: E.pseudo_torsion_polynomial(3) 
     1330            sage: E.division_polynomial_0(3) 
    13571331            3*x^4 - 6*x^2 + 3*x - 1 
    1358             sage: E.pseudo_torsion_polynomial(4) 
     1332            sage: E.division_polynomial_0(4) 
    13591333            2*x^6 - 10*x^4 + 10*x^3 - 10*x^2 + 2*x + 1 
    1360             sage: E.pseudo_torsion_polynomial(5) 
     1334            sage: E.division_polynomial_0(5) 
    13611335            5*x^12 - 62*x^10 + 95*x^9 - 105*x^8 - 60*x^7 + 285*x^6 - 174*x^5 - 5*x^4 - 5*x^3 + 35*x^2 - 15*x + 2 
    1362             sage: E.pseudo_torsion_polynomial(6) 
     1336            sage: E.division_polynomial_0(6) 
    13631337            3*x^16 - 72*x^14 + 168*x^13 - 364*x^12 + 1120*x^10 - 1144*x^9 + 300*x^8 - 540*x^7 + 1120*x^6 - 588*x^5 - 133*x^4 + 252*x^3 - 114*x^2 + 22*x - 1 
    1364             sage: E.pseudo_torsion_polynomial(7) 
     1338            sage: E.division_polynomial_0(7) 
    13651339            7*x^24 - 308*x^22 + 986*x^21 - 2954*x^20 + 28*x^19 + 17171*x^18 - 23142*x^17 + 511*x^16 - 5012*x^15 + 43804*x^14 - 7140*x^13 - 96950*x^12 + 111356*x^11 - 19516*x^10 - 49707*x^9 + 40054*x^8 - 124*x^7 - 18382*x^6 + 13342*x^5 - 4816*x^4 + 1099*x^3 - 210*x^2 + 35*x - 3 
    1366             sage: E.pseudo_torsion_polynomial(8) 
     1340            sage: E.division_polynomial_0(8) 
    13671341            4*x^30 - 292*x^28 + 1252*x^27 - 5436*x^26 + 2340*x^25 + 39834*x^24 - 79560*x^23 + 51432*x^22 - 142896*x^21 + 451596*x^20 - 212040*x^19 - 1005316*x^18 + 1726416*x^17 - 671160*x^16 - 954924*x^15 + 1119552*x^14 + 313308*x^13 - 1502818*x^12 + 1189908*x^11 - 160152*x^10 - 399176*x^9 + 386142*x^8 - 220128*x^7 + 99558*x^6 - 33528*x^5 + 6042*x^4 + 310*x^3 - 406*x^2 + 78*x - 5 
    13681342 
    1369             sage: E.pseudo_torsion_polynomial(18) % E.pseudo_torsion_polynomial(6) == 0 
     1343            sage: E.division_polynomial_0(18) % E.division_polynomial_0(6) == 0 
    13701344            True 
    13711345 
    13721346           An example to illustrate the relationship with torsion points. 
    13731347            sage: F = GF(11) 
    13741348            sage: E = EllipticCurve(F, [0, 2]); E 
    13751349            Elliptic Curve defined by y^2  = x^3 + 2 over Finite Field of size 11 
    1376             sage: f = E.pseudo_torsion_polynomial(5); f 
     1350            sage: f = E.division_polynomial_0(5); f 
    13771351            5*x^12 + x^9 + 8*x^6 + 4*x^3 + 7 
    13781352            sage: f.factor() 
    13791353            (5) * (x^2 + 5) * (x^2 + 2*x + 5) * (x^2 + 5*x + 7) * (x^2 + 7*x + 7) * (x^2 + 9*x + 5) * (x^2 + 10*x + 7) 
    …… 
    13841358           
    13851359            sage: K = GF(11^4, 'a') 
    13861360            sage: X = E.change_ring(K) 
    1387             sage: f = X.pseudo_torsion_polynomial(5) 
     1361            sage: f = X.division_polynomial_0(5) 
    13881362            sage: x_coords = f.roots(multiplicities=False); x_coords 
    13891363            [10*a^3 + 4*a^2 + 5*a + 6, 
    13901364             9*a^3 + 8*a^2 + 10*a + 8, 
    …… 
    14171391          (1 : -1 : 1), 
    14181392          (2 : -5 : 1), 
    14191393          (9 : -33 : 1)] 
    1420           sage: pol=E.pseudo_torsion_polynomial(6) 
     1394          sage: pol=E.division_polynomial_0(6) 
    14211395          sage: xlist=pol.roots(multiplicities=False); xlist 
    14221396          [9, 2, -1/3, -5] 
    14231397          sage: [E.lift_x(x, all=True) for x in xlist] 
    …… 
    14471421             if n == -1: 
    14481422                 answer = 4*x**3 + b2*x**2 + 2*b4*x + b6 
    14491423             elif n == -2: 
    1450                  answer = self.pseudo_torsion_polynomial(-1, x, cache) ** 2 
     1424                 answer = self.division_polynomial_0(-1, x, cache) ** 2 
    14511425             elif n == 1 or n == 2: 
    14521426                 answer = x.parent()(1) 
    14531427             elif n == 3: 
    14541428                 answer = 3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8 
    14551429             elif n == 4: 
    1456                  answer = -self.pseudo_torsion_polynomial(-2, x, cache) + \ 
     1430                 answer = -self.division_polynomial_0(-2, x, cache) + \ 
    14571431                          (6*x**2 + b2*x + b4) * \ 
    1458                           self.pseudo_torsion_polynomial(3, x, cache) 
     1432                          self.division_polynomial_0(3, x, cache) 
    14591433             else: 
    14601434                 raise ValueError, "n must be a positive integer (or -1 or -2)" 
    14611435         else: 
    14621436             if n % 2 == 0: 
    14631437                 m = (n-2) // 2 
    1464                  g_mplus3 = self.pseudo_torsion_polynomial(m+3, x, cache) 
    1465                  g_mplus2 = self.pseudo_torsion_polynomial(m+2, x, cache) 
    1466                  g_mplus1 = self.pseudo_torsion_polynomial(m+1, x, cache) 
    1467                  g_m      = self.pseudo_torsion_polynomial(m,   x, cache) 
    1468                  g_mless1 = self.pseudo_torsion_polynomial(m-1, x, cache) 
     1438                 g_mplus3 = self.division_polynomial_0(m+3, x, cache) 
     1439                 g_mplus2 = self.division_polynomial_0(m+2, x, cache) 
     1440                 g_mplus1 = self.division_polynomial_0(m+1, x, cache) 
     1441                 g_m      = self.division_polynomial_0(m,   x, cache) 
     1442                 g_mless1 = self.division_polynomial_0(m-1, x, cache) 
    14691443                 answer = g_mplus1 * \ 
    14701444                          (g_mplus3 * g_m**2 - g_mless1 * g_mplus2**2) 
    14711445             else: 
    14721446                 m = (n-1) // 2 
    1473                  g_mplus2 = self.pseudo_torsion_polynomial(m+2, x, cache) 
    1474                  g_mplus1 = self.pseudo_torsion_polynomial(m+1, x, cache) 
    1475                  g_m      = self.pseudo_torsion_polynomial(m,   x, cache) 
    1476                  g_mless1 = self.pseudo_torsion_polynomial(m-1, x, cache) 
    1477                  B6_sqr   = self.pseudo_torsion_polynomial(-2, x, cache) 
     1447                 g_mplus2 = self.division_polynomial_0(m+2, x, cache) 
     1448                 g_mplus1 = self.division_polynomial_0(m+1, x, cache) 
     1449                 g_m      = self.division_polynomial_0(m,   x, cache) 
     1450                 g_mless1 = self.division_polynomial_0(m-1, x, cache) 
     1451                 B6_sqr   = self.division_polynomial_0(-2, x, cache) 
    14781452                 if m % 2 == 0: 
    14791453                     answer = B6_sqr * g_mplus2 * g_m**3 - \ 
    14801454                              g_mless1 * g_mplus1**3 
    …… 
    14851459         cache[n] = answer 
    14861460         return answer 
    14871461 
     1462    def two_division_polynomial(self, x = None): 
     1463        r""" 
     1464        Returns the 2-division polynomial of this elliptic curve evaluated at x. 
     1465 
     1466        INPUT: 
     1467             x -- optional ring element to use as the "x" variable. If x 
     1468                  is None, then a new polynomial ring will be constructed over 
     1469                  the base ring of the elliptic curve, and its generator will 
     1470                  be used as x. Note that x does not need to be a generator of 
     1471                  a polynomial ring; any ring element is ok. This permits fast 
     1472                  calculation of the torsion polynomial *evaluated* on any 
     1473                  element of a ring. 
     1474 
     1475        EXAMPLES: 
     1476            sage: E=EllipticCurve('5077a1') 
     1477            sage: E.two_division_polynomial() 
     1478            4*x^3 - 28*x + 25 
     1479            sage: E=EllipticCurve(GF(3^2,'a'),[1,1,1,1,1]) 
     1480            sage: E.two_division_polynomial()           
     1481            x^3 + 2*x^2 + 2 
     1482            sage: E.two_division_polynomial().roots()   
     1483            [(2, 1), (2*a, 1), (a + 2, 1)] 
     1484        """ 
     1485        return self.division_polynomial_0(-1,x) 
     1486    
     1487    def division_polynomial(self, m, x=None, two_torsion_multiplicity=2): 
     1488        r""" 
     1489        Returns the $m$-th division polynomial of this elliptic curve. 
     1490 
     1491        INPUT: 
     1492             m -- positive integer. 
     1493             x -- optional ring element to use as the "x" variable. If x 
     1494                  is None, then a new polynomial ring will be constructed over 
     1495                  the base ring of the elliptic curve, and its generator will 
     1496                  be used as x. Note that x does not need to be a generator of 
     1497                  a polynomial ring; any ring element is ok. This permits fast 
     1498                  calculation of the torsion polynomial *evaluated* on any 
     1499                  element of a ring. 
     1500             two_torsion_multiplicity -- 0,1 or 2 
     1501 
     1502                  If 0: for even $m$ when x is None, a univariate 
     1503                  polynomial over the base ring of the curve is 
     1504                  returned, which omits factors whose roots are the 
     1505                  $x$-coordinates of the $2$-torsion points. 
     1506                  Similarly when $x$ is not none, the evaluation of 
     1507                  such a polynomial at $x$ is returned. 
     1508 
     1509                  If 2: for even $m$ when x is None, a univariate 
     1510                  polynomial over the base ring of the curve is 
     1511                  returned, which includes a factor of degree 3 whose 
     1512                  roots are the $x$-coordinates of the $2$-torsion 
     1513                  points.  Similarly when $x$ is not none, the 
     1514                  evaluation of such a polynomial at $x$ is returned. 
     1515 
     1516                  If 1: when x is None, a bivariate polynomial over 
     1517                  the base ring of the curve is returned, which 
     1518                  includes a factor 2*y+a1*x+a3 which has simple zeros 
     1519                  at the $2$-torsion points.  When $x$ is not none, it 
     1520                  should be a tuple of length 2, and the evaluation of 
     1521                  such a polynomial at $x$ is returned. 
     1522 
     1523        EXAMPLES: 
     1524            sage: E = EllipticCurve([0,0,1,-1,0]) 
     1525            sage: E.division_polynomial(1) 
     1526            1 
     1527            sage: E.division_polynomial(2, two_torsion_multiplicity=0) 
     1528            1 
     1529            sage: E.division_polynomial(2, two_torsion_multiplicity=1) 
     1530            2*y + 1 
     1531            sage: E.division_polynomial(2, two_torsion_multiplicity=2) 
     1532            4*x^3 - 4*x + 1 
     1533            sage: E.division_polynomial(2) 
     1534            4*x^3 - 4*x + 1 
     1535            sage: [E.division_polynomial(3, two_torsion_multiplicity=i) for i in range(3)] 
     1536            [3*x^4 - 6*x^2 + 3*x - 1, 3*x^4 - 6*x^2 + 3*x - 1, 3*x^4 - 6*x^2 + 3*x - 1] 
     1537            sage: [type(E.division_polynomial(3, two_torsion_multiplicity=i)) for i in range(3)] 
     1538            [<class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense'>, 
     1539            <type 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular'>, 
     1540            <class 'sage.rings.polynomial.polynomial_element_generic.Polynomial_rational_dense'>] 
     1541 
     1542            sage: E = EllipticCurve([0, -1, 1, -10, -20])  
     1543            sage: R.<z>=PolynomialRing(QQ)                
     1544            sage: E.division_polynomial(4,z,0) 
     1545            2*z^6 - 4*z^5 - 100*z^4 - 790*z^3 - 210*z^2 - 1496*z - 5821 
     1546            sage: E.division_polynomial(4,z)                            
     1547            8*z^9 - 24*z^8 - 464*z^7 - 2758*z^6 + 6636*z^5 + 34356*z^4 + 53510*z^3 + 99714*z^2 + 351024*z + 459859 
     1548 
     1549        This does not work, since when two_torsion_multiplicity is 1, 
     1550        we compute a bivariate polynomial, and must evaluate at 2 
     1551        tuple of length 2:         
     1552            sage: E.division_polynomial(4,z,1) 
     1553            Traceback (most recent call last): 
     1554            ... 
     1555            ValueError: x should be a tuple of length 2 (or None) when two_torsion_multiplicity is 1 
     1556            sage: R.<z,w>=PolynomialRing(QQ,2)                
     1557            sage: E.division_polynomial(4,(z,w),1).factor() 
     1558            (2*w + 1) * (2*z^6 - 4*z^5 - 100*z^4 - 790*z^3 - 210*z^2 - 1496*z - 5821)        
     1559 
     1560        We can also evaluate this bivariate polynomial at a point: 
     1561            sage: P = E(5,5) 
     1562            sage: E.division_polynomial(4,P,two_torsion_multiplicity=1) 
     1563            -1771561 
     1564        """ 
     1565 
     1566        if not two_torsion_multiplicity in [0,1,2]: 
     1567            raise ValueError, "two_torsion_multiplicity must be 0,1 or 2" 
     1568 
     1569        # Coerce the input m to be an integer 
     1570        m = rings.Integer(m) 
     1571 
     1572        if two_torsion_multiplicity == 0:            
     1573            try: 
     1574                return self.__divpoly0[(m,x)] 
     1575            except AttributeError: 
     1576                self.__divpoly0 = {} 
     1577            except KeyError: 
     1578                pass 
     1579            f = self.division_polynomial_0(m,x) 
     1580            self.__divpoly0[(m,x)] = f 
     1581            return f 
     1582 
     1583        if two_torsion_multiplicity == 1:            
     1584            try: 
     1585                return self.__divpoly1[(m,x)] 
     1586            except AttributeError: 
     1587                self.__divpoly1 = {} 
     1588            except KeyError: 
     1589                pass 
     1590            xy = x 
     1591            R, (x,y) = PolynomialRing(self.base_ring(), 2, 'x,y').objgens() 
     1592            a1,a2,a3,a4,a6 = self.a_invariants() 
     1593            f = self.division_polynomial_0(m,x) 
     1594            if m % 2 == 0: 
     1595                f *= (2*y+a1*x+a3) 
     1596            if xy is None: 
     1597                self.__divpoly1[(m,(x,y))] = f 
     1598                return f 
     1599            else: 
     1600                if isinstance(xy,tuple) and len(xy)==2 or isinstance(xy, ell_point.EllipticCurvePoint_field): 
     1601                    fxy = f(xy[0],xy[1]) 
     1602                    self.__divpoly1[(m,xy)] = fxy 
     1603                    return fxy 
     1604                else: 
     1605                    raise ValueError, "x should be a tuple of length 2 (or None) when two_torsion_multiplicity is 1"                        
     1606 
     1607        if two_torsion_multiplicity == 2:            
     1608            try: 
     1609                return self.__divpoly2[(m,x)] 
     1610            except AttributeError: 
     1611                self.__divpoly2 = {} 
     1612            except KeyError: 
     1613                pass 
     1614            f = self.division_polynomial_0(m,x) 
     1615            if m%2 == 0: 
     1616                f *= self.division_polynomial_0(-1,x)            
     1617            self.__divpoly2[(m,x)] = f 
     1618            return f 
     1619 
     1620    torsion_polynomial = division_polynomial 
     1621 
     1622#     def torsion_polynomial_old(self, n, var='x', i=0): 
     1623#         r""" 
     1624#         Returns the n-th torsion polynomial (a.k.a., division polynomial). 
     1625 
     1626#         INPUT: 
     1627#             n -- non-negative integer 
     1628#             var -- string; the indeterminate of the polynomial (default: x) 
     1629#             i -- integer, either 0 (default) or 1. 
     1630             
     1631#         OUTPUT: 
     1632#             Polynomial -- n-th torsion polynomial, which is a polynomial over 
     1633#                           the base field of the elliptic curve. 
     1634 
     1635#         SEE ALSO: _division_polynomial, full_division_polynomial 
     1636 
     1637#         EXAMPLES: 
     1638#             sage: E = EllipticCurve([0,0,1,-1,0]) 
     1639#             sage: E.torsion_polynomial_old(1) 
     1640#             1 
     1641#             sage: E.torsion_polynomial_old(2) 
     1642#             4*x^3 - 4*x + 1 
     1643#             sage: E.torsion_polynomial_old(3, 'z') 
     1644#             3*z^4 - 6*z^2 + 3*z - 1 
     1645 
     1646#             sage: E = EllipticCurve([0, -1, 1, -10, -20]) 
     1647#             sage: E.torsion_polynomial(0) 
     1648#             0 
     1649#             sage: E.torsion_polynomial(1) 
     1650#             1 
     1651#             sage: E.torsion_polynomial(2) 
     1652#             4*x^3 - 4*x^2 - 40*x - 79 
     1653#             sage: E.torsion_polynomial(3) 
     1654#             3*x^4 - 4*x^3 - 60*x^2 - 237*x - 21 
     1655#             sage: E.torsion_polynomial(4) 
     1656#             8*x^9 - 24*x^8 - 464*x^7 - 2758*x^6 + 6636*x^5 + 34356*x^4 + 53510*x^3 + 99714*x^2 + 351024*x + 459859 
     1657 
     1658#             sage: E = EllipticCurve([-4,0]) 
     1659#             sage: E.torsion_polynomial(2) 
     1660#             4*x^3 - 16*x 
     1661#             sage: E.torsion_polynomial(5) 
     1662#             5*x^12 - 248*x^10 - 1680*x^8 + 19200*x^6 - 32000*x^4 + 51200*x^2 + 4096 
     1663#             sage: E.torsion_polynomial(6) 
     1664#             12*x^19 - 1200*x^17 - 18688*x^15 + 422912*x^13 - 2283520*x^11 + 9134080*x^9 - 27066368*x^7 + 19136512*x^5 + 19660800*x^3 - 3145728*x 
     1665 
     1666#             Check for consistency with division_polynomial_0: 
     1667#             sage: assert all([E.torsion_polynomial(2*m-1)==E.division_polynomial_0(2*m-1) for m in range(1,20)]) 
     1668#             sage: assert all([E.torsion_polynomial(2*m)==E.division_polynomial_0(2*m)*E.division_polynomial_0(-1) for m in range(1,20)]) 
     1669       
     1670#         AUTHOR: David Kohel (kohel@maths.usyd.edu.au), 2005-04-25 
     1671#         """ 
     1672#         n = int(n) 
     1673#         try: 
     1674#             return self.__torsion_polynomial[n] 
     1675#         except AttributeError: 
     1676#             self.__torsion_polynomial = {} 
     1677#         except KeyError: 
     1678#             pass 
     1679#         E = self; i=int(i) 
     1680#         if n < 0: 
     1681#             raise ValueError, "n must be a non-negative integer." 
     1682#         if i > 1 :  
     1683#             raise ValueError, "i must be 0 or 1." 
     1684         
     1685#         R = rings.PolynomialRing(E.base_ring(), var) 
     1686 
     1687# # This is just an abbreviation to make the following code clearer: 
     1688#         tp = lambda m: E.torsion_polynomial(m) 
     1689 
     1690#         if i == 1: 
     1691#             if n == 0: 
     1692#                 f = tp(1) 
     1693#                 E.__torsion_polynomial[n] = f  
     1694#                 return f 
     1695#             else:  
     1696#                 x = R.gen() 
     1697#                 psi2 = tp(2) 
     1698#                 if n%2 == 0:  
     1699#                     f = x * psi2 * (tp(n)//psi2)**2 - tp(n+1) * tp(n-1) 
     1700#                     E.__torsion_polynomial[n] = f 
     1701#                     return f 
     1702#                 else: 
     1703#                     f = x * tp(n)**2 - (tp(n+1)//psi2) * tp(n-1) 
     1704#                     E.__torsion_polynomial[n] = f 
     1705#                     return f 
     1706                 
     1707#         else: 
     1708             
     1709#             if n == 0: return R(0) 
     1710#             if n == 1: return R(1) 
     1711#             x = R.gen() 
     1712#             b2, b4, b6, b8 = E.b_invariants() 
     1713#             psi2 = 4*x**3 + b2*x**2 + 2*b4*x + b6 
     1714#             if n == 2:  
     1715#                 f = psi2 
     1716#                 E.__torsion_polynomial[n] = f; return f 
     1717#             if n == 3: 
     1718#                 f = 3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8 
     1719#                 E.__torsion_polynomial[n] = f; return f 
     1720#             if n == 4: 
     1721#                 f = psi2 * (2*x**6 + b2*x**5 + 5*b4*x**4 + 10*b6*x**3 \ 
     1722#                     + 10*b8*x**2  + (b2*b8 - b4*b6)*x + b4*b8 - b6**2) 
     1723#                 E.__torsion_polynomial[n] = f; return f 
     1724#             if n%2 == 0: 
     1725#                 m = n//2  
     1726#                 if m%2 == 0: 
     1727#                     f = tp(m) * \ 
     1728#                         ((tp(m+2)//psi2) * tp(m-1)**2 - \ 
     1729#                          (tp(m-2)//psi2) * tp(m+1)**2) 
     1730#                     E.__torsion_polynomial[n] = f; return f 
     1731#                 else: 
     1732#                     f = psi2 * tp(m)*( \ 
     1733#                         tp(m+2) * (tp(m-1)//psi2)**2 - \ 
     1734#                         tp(m-2) * (tp(m+1)//psi2)**2) 
     1735#                     E.__torsion_polynomial[n] = f; return f                     
     1736#             else:   
     1737#                 m = n//2  
     1738#                 if m%2 == 0: 
     1739#                     f = psi2 * \ 
     1740#                         tp(m+2) * (tp(m)//psi2)**3 - \ 
     1741#                         tp(m-1) * tp(m+1)**3 
     1742#                     E.__torsion_polynomial[n] = f; return f                     
     1743#                 else: 
     1744#                     f = tp(m+2) * tp(m)**3 - psi2 * \ 
     1745#                         tp(m-1)*(tp(m+1)//psi2)**3 
     1746#                     E.__torsion_polynomial[n] = f; return f 
     1747 
     1748#     def full_division_polynomial(self, m): 
     1749#         """ 
     1750#         Return the $m$-th bivariate division polynomial in $x$ and 
     1751#         $y$.  When $m$ is odd this is exactly the same as the usual 
     1752#         $m$th division polynomial. 
     1753 
     1754#         For the usual division polynomial only in $x$, see the 
     1755#         functions division_polynomial() and 
     1756#         division_polynomial_0(). 
     1757 
     1758#         INPUT: 
     1759#             self -- elliptic curve 
     1760#             m    -- a positive integer 
     1761#         OUTPUT: 
     1762#             a polynomial in two variables $x$, $y$. 
     1763 
     1764#         NOTE: The result is cached.             
     1765 
     1766#         EXAMPLES: 
     1767#         We create a curve and compute the first two full division 
     1768#         polynomials. 
     1769#             sage: E = EllipticCurve([2,3]) 
     1770#             sage: E.full_division_polynomial(1) 
     1771#             1 
     1772#             sage: E.full_division_polynomial(2) 
     1773#             2*y 
     1774 
     1775#         Note that for odd input the full division polynomial is 
     1776#         just the usual division polynomial, but not for even input: 
     1777#             sage: E.division_polynomial(2) 
     1778#             4*x^3 + 8*x + 12 
     1779#             sage: E.full_division_polynomial(3) 
     1780#             3*x^4 + 12*x^2 + 36*x - 4 
     1781#             sage: E.division_polynomial(3) 
     1782#             3*x^4 + 12*x^2 + 36*x - 4 
     1783#             sage: E.full_division_polynomial(4) 
     1784#             4*x^6*y + 40*x^4*y + 240*x^3*y - 80*x^2*y - 96*x*y - 320*y 
     1785#             sage: E.full_division_polynomial(5) 
     1786#             5*x^12 + 124*x^10 + 1140*x^9 - 420*x^8 + 1440*x^7 - 4560*x^6 - 8352*x^5 - 36560*x^4 - 45120*x^3 - 10240*x^2 - 39360*x - 22976 
     1787 
     1788#         TESTS: 
     1789#         We test that the full division polynomial as computed using 
     1790#         the recurrence agrees with the normal division polynomial for 
     1791#         a certain curve and all odd $n$ up to $23$: 
     1792         
     1793#             sage: E = EllipticCurve([23,-105]) 
     1794#             sage: for n in [1,3,..,23]: 
     1795#             ...       assert E.full_division_polynomial(n) == E.division_polynomial(n)         
     1796#         """ 
     1797#         # Coerce the input m to be an integer 
     1798#         m = rings.Integer(m) 
     1799 
     1800#         # Check whether the corresponding poly is cached already 
     1801#         try: 
     1802#             return self.__divpoly2[m] 
     1803#         except AttributeError: 
     1804#             self.__divpoly2 = {} 
     1805#         except KeyError: 
     1806#             pass 
     1807 
     1808#         # Get the invariants of the curve 
     1809#         a1,a2,a3,a4,a6 = self.a_invariants() 
     1810 
     1811#         # Define the polynomial ring that will contain the full 
     1812#         # division polynomial. 
     1813         
     1814#         R, (x,y) = PolynomialRing(self.base_ring(), 2, 'x,y').objgens() 
     1815 
     1816#         # The function division_polynomial_0() gives the correct 
     1817#         # result when m is odd, and all we do in this case is evaluate 
     1818#         # this at the variable x in our bivariate polynomial ring. 
     1819 
     1820#         # For even m, we must multiply the result of 
     1821#         # division_polynomial_0() by the full 2-torsion polynomial 
     1822#         # which is 2*y+a1*x+a3 
     1823 
     1824#         f = self.division_polynomial_0(m,x) 
     1825#         if m % 2 == 0: 
     1826#             f *= (2*y+a1*x+a3) 
     1827 
     1828#         # Cache the result and return it.  
     1829#         self.__divpoly2[m] = f 
     1830#         return f 
    14881831 
    14891832    def multiple_x_numerator(self, n, x=None, cache=None): 
    14901833         r""" 
    14911834         Returns the numerator of the x-coordinate of the nth multiple of 
    14921835         a point, using torsion polynomials (division polynomials). 
    14931836 
    1494          The inputs n, x, cache are as described in pseudo_torsion_polynomial(). 
     1837         The inputs n, x, cache are as described in division_polynomial_0(). 
    14951838 
    14961839         The result is adjusted to be correct for both even and odd n. 
    14971840 
    1498          WARNING: 
    1499            -- There may of course be cancellation between the numerator and 
    1500            the denominator (multiple_x_denominator()). Be careful. For more 
    1501            information on how to avoid cancellation, see Christopher Wuthrich's 
    1502            thesis. 
     1841         WARNING: -- There may of course be cancellation between the 
     1842         numerator and the denominator (multiple_x_denominator()). Be 
     1843         careful. For more information on how to avoid cancellation, 
     1844         see Chris Wuthrich's thesis. 
    15031845 
    15041846         SEE ALSO: 
    15051847           -- multiple_x_denominator() 
    …… 
    15491891         if n < 2: 
    15501892             print "n must be at least 2" 
    15511893 
    1552          self.pseudo_torsion_polynomial( -2, x, cache) 
    1553          self.pseudo_torsion_polynomial(n-1, x, cache) 
    1554          self.pseudo_torsion_polynomial(n  , x, cache) 
    1555          self.pseudo_torsion_polynomial(n+1, x, cache) 
     1894         self.division_polynomial_0( -2, x, cache) 
     1895         self.division_polynomial_0(n-1, x, cache) 
     1896         self.division_polynomial_0(n  , x, cache) 
     1897         self.division_polynomial_0(n+1, x, cache) 
    15561898 
    15571899         if n % 2 == 0: 
    15581900             return x * cache[-1] * cache[n]**2 - cache[n-1] * cache[n+1] 
    …… 
    15651907         Returns the denominator of the x-coordinate of the nth multiple of 
    15661908         a point, using torsion polynomials (division polynomials). 
    15671909 
    1568          The inputs n, x, cache are as described in pseudo_torsion_polynomial(). 
     1910         The inputs n, x, cache are as described in division_polynomial_0(). 
    15691911 
    15701912         The result is adjusted to be correct for both even and odd n. 
    15711913 
    …… 
    15931935         if n < 2: 
    15941936             print "n must be at least 2" 
    15951937 
    1596          self.pseudo_torsion_polynomial(-2, x, cache) 
    1597          self.pseudo_torsion_polynomial(n , x, cache) 
     1938         self.division_polynomial_0(-2, x, cache) 
     1939         self.division_polynomial_0(n , x, cache) 
    15981940 
    15991941         if n % 2 == 0: 
    16001942             return cache[-1] * cache[n]**2 
    16011943         else: 
    16021944             return cache[n]**2 
    16031945 
    1604  
    1605     def torsion_polynomial(self, n, var='x', i=0): 
    1606         r""" 
    1607         Returns the n-th torsion polynomial (a.k.a., division polynomial). 
    1608  
    1609         INPUT: 
    1610             n -- non-negative integer 
    1611             var -- string; the indeterminate of the polynomial (default: x) 
    1612             i -- integer, either 0 (default) or 1. 
    1613              
    1614         OUTPUT: 
    1615             Polynomial -- n-th torsion polynomial, which is a polynomial over 
    1616                           the base field of the elliptic curve. 
    1617  
    1618         SEE ALSO: pseudo_torsion_polynomial, full_division_polynomial 
    1619  
    1620         ALIASES: division_polynomial, torsion_polynomial 
    1621  
    1622         EXAMPLES: 
    1623             sage: E = EllipticCurve([0,0,1,-1,0]) 
    1624             sage: E.division_polynomial(1) 
    1625             1 
    1626             sage: E.division_polynomial(2) 
    1627             4*x^3 - 4*x + 1 
    1628             sage: E.division_polynomial(3, 'z') 
    1629             3*z^4 - 6*z^2 + 3*z - 1 
    1630  
    1631             sage: E = EllipticCurve([0, -1, 1, -10, -20]) 
    1632             sage: E.torsion_polynomial(0) 
    1633             0 
    1634             sage: E.torsion_polynomial(1) 
    1635             1 
    1636             sage: E.torsion_polynomial(2) 
    1637             4*x^3 - 4*x^2 - 40*x - 79 
    1638             sage: E.torsion_polynomial(3) 
    1639             3*x^4 - 4*x^3 - 60*x^2 - 237*x - 21 
    1640             sage: E.torsion_polynomial(4) 
    1641             8*x^9 - 24*x^8 - 464*x^7 - 2758*x^6 + 6636*x^5 + 34356*x^4 + 53510*x^3 + 99714*x^2 + 351024*x + 459859 
    1642  
    1643             sage: E = EllipticCurve([-4,0]) 
    1644             sage: E.torsion_polynomial(2) 
    1645             4*x^3 - 16*x 
    1646             sage: E.torsion_polynomial(5) 
    1647             5*x^12 - 248*x^10 - 1680*x^8 + 19200*x^6 - 32000*x^4 + 51200*x^2 + 4096 
    1648             sage: E.torsion_polynomial(6) 
    1649             12*x^19 - 1200*x^17 - 18688*x^15 + 422912*x^13 - 2283520*x^11 + 9134080*x^9 - 27066368*x^7 + 19136512*x^5 + 19660800*x^3 - 3145728*x 
    1650  
    1651             Check for consistency with pseudo_torsion_polynomial: 
    1652             sage: assert all([E.torsion_polynomial(2*m-1)==E.pseudo_torsion_polynomial(2*m-1) for m in range(1,20)]) 
    1653             sage: assert all([E.torsion_polynomial(2*m)==E.pseudo_torsion_polynomial(2*m)*E.pseudo_torsion_polynomial(-1) for m in range(1,20)]) 
    1654        
    1655         AUTHOR: David Kohel (kohel@maths.usyd.edu.au), 2005-04-25 
    1656         """ 
    1657         n = int(n) 
    1658         try: 
    1659             return self.__torsion_polynomial[n] 
    1660         except AttributeError: 
    1661             self.__torsion_polynomial = {} 
    1662         except KeyError: 
    1663             pass 
    1664         E = self; i=int(i) 
    1665         if n < 0: 
    1666             raise ValueError, "n must be a non-negative integer." 
    1667         if i > 1 :  
    1668             raise ValueError, "i must be 0 or 1." 
    1669          
    1670         R = rings.PolynomialRing(E.base_ring(), var) 
    1671  
    1672 # This is just an abbreviation to make the following code clearer: 
    1673         tp = lambda m: E.torsion_polynomial(m) 
    1674  
    1675         if i == 1: 
    1676             if n == 0: 
    1677                 f = tp(1) 
    1678                 E.__torsion_polynomial[n] = f  
    1679                 return f 
    1680             else:  
    1681                 x = R.gen() 
    1682                 psi2 = tp(2) 
    1683                 if n%2 == 0:  
    1684                     f = x * psi2 * (tp(n)//psi2)**2 - tp(n+1) * tp(n-1) 
    1685                     E.__torsion_polynomial[n] = f 
    1686                     return f 
    1687                 else: 
    1688                     f = x * tp(n)**2 - (tp(n+1)//psi2) * tp(n-1) 
    1689                     E.__torsion_polynomial[n] = f 
    1690                     return f 
    1691                  
    1692         else: 
    1693              
    1694             if n == 0: return R(0) 
    1695             if n == 1: return R(1) 
    1696             x = R.gen() 
    1697             b2, b4, b6, b8 = E.b_invariants() 
    1698             psi2 = 4*x**3 + b2*x**2 + 2*b4*x + b6 
    1699             if n == 2:  
    1700                 f = psi2 
    1701                 E.__torsion_polynomial[n] = f; return f 
    1702             if n == 3: 
    1703                 f = 3*x**4 + b2*x**3 + 3*b4*x**2 + 3*b6*x + b8 
    1704                 E.__torsion_polynomial[n] = f; return f 
    1705             if n == 4: 
    1706                 f = psi2 * (2*x**6 + b2*x**5 + 5*b4*x**4 + 10*b6*x**3 \ 
    1707                     + 10*b8*x**2  + (b2*b8 - b4*b6)*x + b4*b8 - b6**2) 
    1708                 E.__torsion_polynomial[n] = f; return f 
    1709             if n%2 == 0: 
    1710                 m = n//2  
    1711                 if m%2 == 0: 
    1712                     f = tp(m) * \ 
    1713                         ((tp(m+2)//psi2) * tp(m-1)**2 - \ 
    1714                          (tp(m-2)//psi2) * tp(m+1)**2) 
    1715                     E.__torsion_polynomial[n] = f; return f 
    1716                 else: 
    1717                     f = psi2 * tp(m)*( \ 
    1718                         tp(m+2) * (tp(m-1)//psi2)**2 - \ 
    1719                         tp(m-2) * (tp(m+1)//psi2)**2) 
    1720                     E.__torsion_polynomial[n] = f; return f                     
    1721             else:   
    1722                 m = n//2  
    1723                 if m%2 == 0: 
    1724                     f = psi2 * \ 
    1725                         tp(m+2) * (tp(m)//psi2)**3 - \ 
    1726                         tp(m-1) * tp(m+1)**3 
    1727                     E.__torsion_polynomial[n] = f; return f                     
    1728                 else: 
    1729                     f = tp(m+2) * tp(m)**3 - psi2 * \ 
    1730                         tp(m-1)*(tp(m+1)//psi2)**3 
    1731                     E.__torsion_polynomial[n] = f; return f 
    1732  
    1733     division_polynomial = torsion_polynomial 
    1734  
    1735     def full_division_polynomial(self, m): 
    1736         """ 
    1737         Return the $m$-th bivariate division polynomial in $x$ and 
    1738         $y$.  When $m$ is odd this is exactly the same as the usual 
    1739         $m$th division polynomial. 
    1740  
    1741         For the usual division polynomial only in $x$, see the 
    1742         functions division_polynomial() and 
    1743         pseudo_torsion_polynomial(). 
    1744  
    1745         INPUT: 
    1746             self -- elliptic curve 
    1747             m    -- a positive integer 
    1748         OUTPUT: 
    1749             a polynomial in two variables $x$, $y$. 
    1750  
    1751         NOTE: The result is cached.             
    1752  
    1753         EXAMPLES: 
    1754         We create a curve and compute the first two full division 
    1755         polynomials. 
    1756             sage: E = EllipticCurve([2,3]) 
    1757             sage: E.full_division_polynomial(1) 
    1758             1 
    1759             sage: E.full_division_polynomial(2) 
    1760             2*y 
    1761  
    1762         Note that for odd input the full division polynomial is 
    1763         just the usual division polynomial, but not for even input: 
    1764             sage: E.division_polynomial(2) 
    1765             4*x^3 + 8*x + 12 
    1766             sage: E.full_division_polynomial(3) 
    1767             3*x^4 + 12*x^2 + 36*x - 4 
    1768             sage: E.division_polynomial(3) 
    1769             3*x^4 + 12*x^2 + 36*x - 4 
    1770             sage: E.full_division_polynomial(4) 
    1771             4*x^6*y + 40*x^4*y + 240*x^3*y - 80*x^2*y - 96*x*y - 320*y 
    1772             sage: E.full_division_polynomial(5) 
    1773             5*x^12 + 124*x^10 + 1140*x^9 - 420*x^8 + 1440*x^7 - 4560*x^6 - 8352*x^5 - 36560*x^4 - 45120*x^3 - 10240*x^2 - 39360*x - 22976 
    1774  
    1775         TESTS: 
    1776         We test that the full division polynomial as computed using 
    1777         the recurrence agrees with the normal division polynomial for 
    1778         a certain curve and all odd $n$ up to $23$: 
    1779          
    1780             sage: E = EllipticCurve([23,-105]) 
    1781             sage: for n in [1,3,..,23]: 
    1782             ...       assert E.full_division_polynomial(n) == E.division_polynomial(n)         
    1783         """ 
    1784         # Coerce the input m to be an integer 
    1785         m = rings.Integer(m) 
    1786  
    1787         # Check whether the corresponding poly is cached already 
    1788         try: 
    1789             return self.__divpoly2[m] 
    1790         except AttributeError: 
    1791             self.__divpoly2 = {} 
    1792         except KeyError: 
    1793             pass 
    1794  
    1795         # Get the invariants of the curve 
    1796