/sage/schemes/elliptic_curves/ell_point.py

Large files files are truncated, but you can click here to view the full file

# -*- coding: utf-8 -*-
r"""
Points on elliptic curves

The base class ``EllipticCurvePoint_field``, derived from
``AdditiveGroupElement``, provides support for points on elliptic
curves defined over general fields.  The derived classes
``EllipticCurvePoint_number_field`` and
``EllipticCurvePoint_finite_field`` provide further support for point
on curves defined over number fields (including the rational field
`\QQ`) and over finite fields.  Although there is no special
class for points over `\QQ`, there is currently greater
functionality implemented over `\QQ` than over other number
fields.

The class ``EllipticCurvePoint``, which is based on
``SchemeMorphism_projective_coordinates_ring``, currently has little
extra functionality.

EXAMPLES:

An example over `\QQ`::

    sage: E = EllipticCurve('389a1')
    sage: P = E(-1,1); P
    (-1 : 1 : 1)
    sage: Q = E(0,-1); Q
    (0 : -1 : 1)
    sage: P+Q
    (4 : 8 : 1)
    sage: P-Q
    (1 : 0 : 1)
    sage: 3*P-5*Q
    (328/361 : -2800/6859 : 1)

An example over a number field::

    sage: K.<i> = QuadraticField(-1)
    sage: E = EllipticCurve(K,[1,0,0,0,-1])
    sage: P = E(0,i); P
    (0 : i : 1)
    sage: P.order()
    +Infinity
    sage: 101*P-100*P==P
    True

An example over a finite field::

    sage: K.<a> = GF(101^3)
    sage: E = EllipticCurve(K,[1,0,0,0,-1])
    sage: P = E(40*a^2 + 69*a + 84 , 58*a^2 + 73*a + 45)
    sage: P.order()
    1032210
    sage: E.cardinality()
    1032210

Arithmetic with a point over an extension of a finite field::

    sage: k.<a> = GF(5^2)
    sage: E = EllipticCurve(k,[1,0]); E
    Elliptic Curve defined by y^2 = x^3 + x over Finite Field in a of size 5^2
    sage: P = E([a,2*a+4])
    sage: 5*P
    (2*a + 3 : 2*a : 1)
    sage: P*5
    (2*a + 3 : 2*a : 1)
    sage: P + P + P + P + P
    (2*a + 3 : 2*a : 1)

::

    sage: F = Zmod(3)
    sage: E = EllipticCurve(F,[1,0]);
    sage: P = E([2,1])
    sage: import sys
    sage: n = sys.maxint
    sage: P*(n+1)-P*n == P
    True

Arithmetic over `\ZZ/N\ZZ` with composite `N` is supported.  When an
operation tries to invert a non-invertible element, a
ZeroDivisionError is raised and a factorization of the modulus appears
in the error message::
    
    sage: N = 1715761513
    sage: E = EllipticCurve(Integers(N),[3,-13])
    sage: P = E(2,1)
    sage: LCM([2..60])*P
    Traceback (most recent call last):
    ...
    ZeroDivisionError: Inverse of 1520944668 does not exist (characteristic = 1715761513 = 26927*63719)


AUTHORS:

- William Stein (2005) -- Initial version
- Robert Bradshaw et al....
- John Cremona (Feb 2008) -- Point counting and group structure for
   non-prime fields, Frobenius endomorphism and order, elliptic logs
- John Cremona (Aug 2008) -- Introduced ``EllipticCurvePoint_number_field`` class 
- Tobias Nagel, Michael Mardaus, John Cremona (Dec 2008) -- `p`-adic elliptic logarithm over `\QQ`
- David Hansen (Jan 2009) -- Added ``weil_pairing`` function to ``EllipticCurvePoint_finite_field`` class
"""

#*****************************************************************************
#       Copyright (C) 2005 William Stein <wstein@gmail.com>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#    This code is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#    General Public License for more details.
#
#  The full text of the GPL is available at:
#
#                  http://www.gnu.org/licenses/
#*****************************************************************************

import math

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
from sage.rings.padics.precision_error import PrecisionError

import ell_generic
import sage.rings.all as rings
import sage.rings.arith as arith
import sage.misc.misc as misc
from sage.groups.all import AbelianGroup
import sage.groups.generic as generic
from sage.libs.pari.all import pari, PariError

from sage.structure.sequence  import Sequence
from sage.schemes.generic.morphism import (SchemeMorphism_projective_coordinates_ring,
                                           SchemeMorphism_abelian_variety_coordinates_field,
                                           is_SchemeMorphism, SchemeMorphism_coordinates)

import sage.schemes.generic.scheme as scheme
from constructor import EllipticCurve
                                           
oo = rings.infinity       # infinity

class EllipticCurvePoint(SchemeMorphism_projective_coordinates_ring):
    """
    A point on an elliptic curve.
    """
    def __cmp__(self, other):
        """
        Standard comparison function for points on elliptic curves, to
        allow sorting and equality testing.

        EXAMPLES:
            sage: E=EllipticCurve(QQ,[1,1])
            sage: P=E(0,1)
            sage: P.order()
            +Infinity
            sage: Q=P+P
            sage: P==Q
            False
            sage: Q+Q == 4*P
            True
        """
        if isinstance(other, (int, long, rings.Integer)) and other == 0:
            if self.is_zero():
                return 0
            else:
                return -1
        return SchemePoint_projective_abelian_scheme.__cmp__(self, other)

class EllipticCurvePoint_field(AdditiveGroupElement): # SchemeMorphism_abelian_variety_coordinates_field):
    """
    A point on an elliptic curve over a field.  The point has coordinates
    in the base field.

    EXAMPLES::

        sage: E = EllipticCurve('37a')
        sage: E([0,0])
        (0 : 0 : 1)
        sage: E(0,0)               # brackets are optional
        (0 : 0 : 1)
        sage: E([GF(5)(0), 0])     # entries are coerced
        (0 : 0 : 1)      

        sage: E(0.000, 0)
        (0 : 0 : 1)
         
        sage: E(1,0,0)  
        Traceback (most recent call last):
        ...
        TypeError: Coordinates [1, 0, 0] do not define a point on
        Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field

    ::

        sage: E = EllipticCurve([0,0,1,-1,0])
        sage: S = E(QQ); S
        Abelian group of points on Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field

        sage: K.<i>=NumberField(x^2+1)
        sage: E=EllipticCurve(K,[0,1,0,-160,308])
        sage: P=E(26,-120)
        sage: Q=E(2+12*i,-36+48*i)
        sage: P.order() == Q.order() == 4  # long time (3s)
        True
        sage: 2*P==2*Q
        False

    ::
        
        sage: K.<t>=FractionField(PolynomialRing(QQ,'t'))
        sage: E=EllipticCurve([0,0,0,0,t^2])    
        sage: P=E(0,t)
        sage: P,2*P,3*P
        ((0 : t : 1), (0 : -t : 1), (0 : 1 : 0))


    TESTS::

        sage: loads(S.dumps()) == S
        True
        sage: E = EllipticCurve('37a')
        sage: P = E(0,0); P
        (0 : 0 : 1)
        sage: loads(P.dumps()) == P
        True
        sage: T = 100*P
        sage: loads(T.dumps()) == T
        True

    Test pickling an elliptic curve that has known points on it::

        sage: e = EllipticCurve([0, 0, 1, -1, 0]); g = e.gens(); loads(dumps(e)) == e
        True
    """
    def __init__(self, curve, v, check=True):
        """
        Constructor for a point on an elliptic curve.

        INPUT:

        - curve -- an elliptic curve
        - v -- data determining a point (another point, the integer
                 0, or a tuple of coordinates)

        EXAMPLE::

            sage: E = EllipticCurve('43a')
            sage: P = E([2, -4, 2]); P
            (1 : -2 : 1)
            sage: P == E([1,-2])
            True
            sage: P = E(0); P
            (0 : 1 : 0)
            sage: P=E(2, -4, 2); P
            (1 : -2 : 1)
        """
        point_homset = curve.point_homset()
        AdditiveGroupElement.__init__(self, point_homset)

        if check:
            # mostly from SchemeMorphism_projective_coordinates_field
            d = point_homset.codomain().ambient_space().ngens()
            if is_SchemeMorphism(v) or isinstance(v, EllipticCurvePoint_field):
                v = list(v)
            elif v == 0:
                self._coords = Sequence((0,1,0), point_homset.value_ring())
                return
            if not isinstance(v,(list,tuple)):
                raise TypeError, \
                      "Argument v (= %s) must be a scheme point, list, or tuple."%str(v)
            if len(v) != d and len(v) != d-1:
                raise TypeError, "v (=%s) must have %s components"%(v, d)
            v = Sequence(v, point_homset.value_ring())
            if len(v) == d-1:     # very common special case
                v.append(v.universe()(1))
            
            n = len(v)
            all_zero = True
            for i in range(n):
                c = v[n-1-i]
                if c:
                    all_zero = False
                    if c == 1:
                        break
                    for j in range(n-i):
                        v[j] /= c
                    break
            if all_zero:
                raise ValueError, "%s does not define a valid point since all entries are 0"%repr(v)

            x, y, z = v
            if z == 0:
                test = x
            else:
                a1, a2, a3, a4, a6 = curve.ainvs()
                test = y**2 + (a1*x+a3)*y - (((x+a2)*x+a4)*x+a6)
            if not test == 0:
                raise TypeError, "Coordinates %s do not define a point on %s"%(list(v),curve)
#            point_homset.codomain()._check_satisfies_equations(v)
                
        self._coords = v
        
        
    def _repr_(self):
        """
        Return a string representation of this point.

        EXAMPLE::

            sage: E = EllipticCurve('39a')
            sage: P = E([-2, 1, 1])
            sage: P._repr_()
            '(-2 : 1 : 1)'
        """
        return self.codomain().ambient_space()._repr_generic_point(self._coords)

    def _latex_(self):
        """
        Return a LaTeX representation of this point.

        EXAMPLE::

            sage: E = EllipticCurve('40a')
            sage: P = E([3, 0])
            sage: P._latex_()
            '\\left(3 : 0 : 1\\right)'
        """
        return self.codomain().ambient_space()._latex_generic_point(self._coords)

    def __getitem__(self, n):
        """
        Return the n'th coordinate of this point.

        EXAMPLE::

            sage: E = EllipticCurve('42a')
            sage: P = E([-17, -51, 17])
            sage: [P[i] for i in [2,1,0]]
            [1, -3, -1]
        """
        return self._coords[n]

    def __iter__(self):
        """
        Return the coordinates of this point as a list.

        EXAMPLE::

            sage: E = EllipticCurve('37a')
            sage: list(E([0,0]))
            [0, 0, 1]
        """
        return iter(self._coords)

    def __tuple__(self):
        """
        Return the coordinates of this point as a tuple.

        EXAMPLE::

            sage: E = EllipticCurve('44a')
            sage: P = E([1, -2, 1])
            sage: P.__tuple__()
            (1, -2, 1)
        """
        return tuple(self._coords) # Warning: _coords is a list!
    
    def __cmp__(self, other):
        """
        Comparison function for points to allow sorting and equality testing.

        EXAMPLES::

            sage: E = EllipticCurve('45a')
            sage: P = E([2, -1, 1])
            sage: P == E(0)
            False
            sage: P+P == E(0)
            True
        """
        if not isinstance(other, EllipticCurvePoint_field):
            try:
                other = self.codomain().ambient_space()(other)
            except TypeError:
                return -1
        return cmp(self._coords, other._coords)

    def _pari_(self):
        r"""
        Converts this point to PARI format.
        
        EXAMPLES::
            
            sage: E = EllipticCurve([0,0,0,3,0])
            sage: O = E(0)
            sage: P = E.point([1,2])
            sage: O._pari_()
            [0]
            sage: P._pari_()
            [1, 2]
        
        The following implicitly calls O._pari_() and P._pari_()::
            
            sage: pari(E).elladd(O,P)
            [1, 2]
        
        TESTS::
            
        Try the same over a finite field::
            
            sage: E = EllipticCurve(GF(11), [0,0,0,3,0])
            sage: O = E(0)
            sage: P = E.point([1,2])
            sage: O._pari_()
            [0]
            sage: P._pari_()
            [Mod(1, 11), Mod(2, 11)]
            sage: pari(E).elladd(O,P)
            [Mod(1, 11), Mod(2, 11)]
        """
        if self[2]:
            return pari([self[0]/self[2], self[1]/self[2]])
        else:
            return pari([0])

    def scheme(self):
        """
        Return the scheme of this point, i.e., the curve it is on.
        This is synonymous with curve() which is perhaps more
        intuitive.
        
        EXAMPLES::

            sage: E=EllipticCurve(QQ,[1,1])
            sage: P=E(0,1)
            sage: P.scheme()
            Elliptic Curve defined by y^2 = x^3 + x + 1 over Rational Field
            sage: P.scheme() == P.curve()
            True
            sage: K.<a>=NumberField(x^2-3,'a')
            sage: P=E.base_extend(K)(1,a)
            sage: P.scheme()
            Elliptic Curve defined by y^2 = x^3 + x + 1 over Number Field in a with defining polynomial x^2 - 3
        """
        #The following text is just not true: it applies to the class
        #EllipticCurvePoint, which appears to be never used, but does
        #not apply to EllipticCurvePoint_field which is simply derived
        #from AdditiveGroupElement.
        #
        #"Technically, points on curves in Sage are scheme maps from
        #  the domain Spec(F) where F is the base field of the curve to
        #  the codomain which is the curve.  See also domain() and
        #  codomain()."

        return self.codomain()

    def domain(self):
        """
        Return the domain of this point, which is `Spec(F)` where `F` is
        the field of definition.

        EXAMPLES::

            sage: E=EllipticCurve(QQ,[1,1])
            sage: P=E(0,1)
            sage: P.domain()
            Spectrum of Rational Field
            sage: K.<a>=NumberField(x^2-3,'a')
            sage: P=E.base_extend(K)(1,a)
            sage: P.domain()
            Spectrum of Number Field in a with defining polynomial x^2 - 3
       """
        return self.parent().domain()
        
    def codomain(self):
        """
        Return the codomain of this point, which is the curve it is
        on.  Synonymous with curve() which is perhaps more intuitive.

        EXAMPLES::

            sage: E=EllipticCurve(QQ,[1,1])
            sage: P=E(0,1)
            sage: P.domain()
            Spectrum of Rational Field
            sage: K.<a>=NumberField(x^2-3,'a')
            sage: P=E.base_extend(K)(1,a)
            sage: P.codomain()
            Elliptic Curve defined by y^2 = x^3 + x + 1 over Number Field in a with defining polynomial x^2 - 3
            sage: P.codomain() == P.curve()
            True
       """
        return self.parent().codomain()

    def order(self):
        r"""
        Return the order of this point on the elliptic curve.

        If the point is zero, returns 1, otherwise raise a
        NotImplementedError.

        For curves over number fields and finite fields, see below.

        .. note::

           :meth:`additive_order` is a synonym for :meth:`order`

        EXAMPLE::

            sage: K.<t>=FractionField(PolynomialRing(QQ,'t'))
            sage: E=EllipticCurve([0,0,0,-t^2,0])   
            sage: P=E(t,0) 
            sage: P.order()                                  
            Traceback (most recent call last):
            ...
            NotImplementedError: Computation of order of a point not implemented over general fields.
            sage: E(0).additive_order()
            1
            sage: E(0).order() == 1
            True

        """
        if self.is_zero():
            return rings.Integer(1)
        raise NotImplementedError, "Computation of order of a point not implemented over general fields."

    additive_order = order
    
    def curve(self):
        """
        Return the curve that this point is on.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1])
            sage: P.curve()
            Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field
        """
        return self.scheme()

    def __nonzero__(self):
        """
        Return True if this is not the zero point on the curve.

        EXAMPLES::

            sage: E = EllipticCurve('37a')
            sage: P = E(0); P
            (0 : 1 : 0)
            sage: P.is_zero()
            True
            sage: P = E.gens()[0]
            sage: P.is_zero()
            False
        """
        return bool(self[2])

    def has_finite_order(self):
        """
        Return True if this point has finite additive order as an element
        of the group of points on this curve.

        For fields other than number fields and finite fields, this is
        NotImplemented unless self.is_zero().

        EXAMPLES::

            sage: K.<t>=FractionField(PolynomialRing(QQ,'t'))
            sage: E=EllipticCurve([0,0,0,-t^2,0])   
            sage: P = E(0)
            sage: P.has_finite_order()
            True
            sage: P=E(t,0) 
            sage: P.has_finite_order()                                  
            Traceback (most recent call last):
            ...
            NotImplementedError: Computation of order of a point not implemented over general fields.
            sage: (2*P).is_zero()
            True
        """
        if self.is_zero(): return True
        return self.order() != oo
    
    is_finite_order = has_finite_order # for backward compatibility

    def has_infinite_order(self):
        """
        Return True if this point has infinite additive order as an element
        of the group of points on this curve.

        For fields other than number fields and finite fields, this is
        NotImplemented unless self.is_zero().

        EXAMPLES::

            sage: K.<t>=FractionField(PolynomialRing(QQ,'t'))
            sage: E=EllipticCurve([0,0,0,-t^2,0])   
            sage: P = E(0)
            sage: P.has_infinite_order()
            False
            sage: P=E(t,0) 
            sage: P.has_infinite_order()                                  
            Traceback (most recent call last):
            ...
            NotImplementedError: Computation of order of a point not implemented over general fields.
            sage: (2*P).is_zero()
            True
        """
        if self.is_zero(): return False
        return self.order() == oo
    
    def plot(self, **args):
        """
        Plot this point on an elliptic curve.

        INPUT:

        - ``**args`` -- all arguments get passed directly onto the point
          plotting function.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1])
            sage: P.plot(pointsize=30, rgbcolor=(1,0,0))
        """
        if self.is_zero():
            return plot.text("$\\infty$", (-3,3), **args)
                   
        else:
            return plot.point((self[0], self[1]), **args)

    def _add_(self, right):
        """
        Add self to right.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1]); Q = E([0,0])
            sage: P + Q
            (1 : 0 : 1)
            sage: P._add_(Q) == P + Q
            True

        Example to show that bug \#4820 is fixed::

            sage: [type(c) for c in 2*EllipticCurve('37a1').gen(0)]
            [<type 'sage.rings.rational.Rational'>,
            <type 'sage.rings.rational.Rational'>,
            <type 'sage.rings.rational.Rational'>]
        """
        # Use Prop 7.1.7 of Cohen "A Course in Computational Algebraic Number Theory"
        if self.is_zero():
            return right
        if right.is_zero():
            return self
        E = self.curve()
        a1, a2, a3, a4, a6 = E.ainvs()
        x1, y1 = self[0], self[1]
        x2, y2 = right[0], right[1]
        if x1 == x2 and y1 == -y2 - a1*x2 - a3:
            return E(0) # point at infinity

        if x1 == x2 and y1 == y2:
            try:
                m = (3*x1*x1 + 2*a2*x1 + a4 - a1*y1) / (2*y1 + a1*x1 + a3)
            except ZeroDivisionError:
                R = E.base_ring()
                if R.is_finite():
                    N = R.characteristic()
                    from sage.rings.all import ZZ
                    N1 = N.gcd(ZZ(2*y1 + a1*x1 + a3))
                    N2 = N//N1
                    raise ZeroDivisionError, "Inverse of %s does not exist (characteristic = %s = %s*%s)"%(2*y1 + a1*x1 + a3, N,N1,N2)
                else:
                    raise ZeroDivisionError, "Inverse of %s does not exist"%(2*y1 + a1*x1 + a3)
        else:
            try:
                m = (y1-y2)/(x1-x2)
            except ZeroDivisionError:
                R = E.base_ring()
                if R.is_finite():
                    N = R.characteristic()
                    from sage.rings.all import ZZ
                    N1 = N.gcd(ZZ(x1-x2))
                    N2 = N//N1
                    raise ZeroDivisionError, "Inverse of %s does not exist (characteristic = %s = %s*%s)"%(x1-x2, N,N1,N2)
                else:
                    raise ZeroDivisionError, "Inverse of %s does not exist"%(x1-x2)

        x3 = -x1 - x2 - a2 + m*(m+a1)
        y3 = -y1 - a3 - a1*x3 + m*(x1-x3)
        # See \#4820 for why we need to coerce 1 into the base ring here:
        return E.point([x3, y3, E.base_ring()(1)], check=False)

    def _sub_(self, right):
        """
        Subtract right from  self.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1]); Q = E([0,0])
            sage: P - Q
            (4 : 8 : 1)
            sage: P - Q == P._sub_(Q)
            True
            sage: (P - Q) + Q
            (-1 : 1 : 1)
            sage: P
            (-1 : 1 : 1)
        """
        return self + (-right)

    def __neg__(self):
        """
        Return the additive inverse of this point.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1])
            sage: Q = -P; Q
            (-1 : -2 : 1)
            sage: Q + P
            (0 : 1 : 0)

        Example to show that bug \#4820 is fixed::

            sage: [type(c) for c in -EllipticCurve('37a1').gen(0)]
            [<type 'sage.rings.rational.Rational'>,
            <type 'sage.rings.rational.Rational'>,
            <type 'sage.rings.rational.Rational'>]
        """
        if self.is_zero():
            return self
        E, x, y = self.curve(), self[0], self[1]
        # See \#4820 for why we need to coerce 1 into the base ring here:
        return E.point([x, -y - E.a1()*x - E.a3(), E.base_ring()(1)], check=False)

    def xy(self):
        """
        Return the `x` and `y` coordinates of this point, as a 2-tuple.
        If this is the point at infinity a ZeroDivisionError is raised.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: P = E([-1,1])
            sage: P.xy()
            (-1, 1)
            sage: Q = E(0); Q
            (0 : 1 : 0)
            sage: Q.xy()
            Traceback (most recent call last):
            ...
            ZeroDivisionError: Rational division by zero
        """
        if self[2] == 1:
            return self[0], self[1]
        else:
            return self[0]/self[2], self[1]/self[2]

    def is_divisible_by(self, m):
        """
        Return True if there exists a point `Q` defined over the same
        field as self such that `mQ` == self.

        INPUT:

        - ``m`` -- a positive integer.

        OUTPUT:

        (bool) -- True if there is a solution, else False.

        .. warning::

           This function usually triggers the computation of the
           `m`-th division polynomial of the associated elliptic
           curve, which will be expensive if `m` is large, though it
           will be cached for subsequent calls with the same `m`.

        EXAMPLES::

            sage: E = EllipticCurve('389a')
            sage: Q = 5*E(0,0); Q
            (-2739/1444 : -77033/54872 : 1)
            sage: Q.is_divisible_by(4)
            False
            sage: Q.is_divisible_by(5)
            True

        A finite field example::

            sage: E = EllipticCurve(GF(101),[23,34])
            sage: E.cardinality().factor()
            2 * 53
            sage: Set([T.order() for T in E.points()])
            {1, 106, 2, 53}
            sage: len([T for T in E.points() if T.is_divisible_by(2)])
            53
            sage: len([T for T in E.points() if T.is_divisible_by(3)])
            106

        TESTS:

        This shows that the bug reported at #10076 is fixed::

            sage: K = QuadraticField(8,'a')
            sage: E = EllipticCurve([K(0),0,0,-1,0])
            sage: P = E([-1,0])
            sage: P.is_divisible_by(2)
            False
            sage: P.division_points(2)
            []

        Note that it is not sufficient to test that
        ``self.division_points(m,poly_only=True)`` has roots::

            sage: P.division_points(2, poly_only=True).roots()
            [(1/2*a - 1, 1), (-1/2*a - 1, 1)]

            sage: tor = E.torsion_points(); len(tor)
            8
            sage: [T.order() for T in tor]
            [2, 4, 4, 2, 4, 1, 4, 2]
            sage: all([T.is_divisible_by(3) for T in tor])
            True
            sage: Set([T for T in tor if T.is_divisible_by(2)])
            {(0 : 1 : 0), (1 : 0 : 1)}
            sage: Set([2*T for T in tor])
            {(0 : 1 : 0), (1 : 0 : 1)}


        """
        # Coerce the input m to an integer
        m = rings.Integer(m)

        # Check for trivial cases of m = 1, -1 and 0.
        if m == 1 or m == -1:
            return True
        if m == 0:
            return self == 0 # then m*self=self for all m!
        m = m.abs()

        # Now the following line would of course be correct, but we
        # work harder to be more efficient:
        # return len(self.division_points(m)) > 0

        P = self

        # If P has finite order n and gcd(m,n)=1 then the result is
        # True.  However, over general fields computing P.order() is
        # not implemented.

        try:
            n = P.order()
            if not n == oo:
                if m.gcd(n)==1:
                    return True
        except NotImplementedError:
            pass

        P_is_2_torsion = (P==-P)
        g = P.division_points(m, poly_only=True)

        if not P_is_2_torsion:
            # In this case deg(g)=m^2, and each root in K lifts to two
            # points Q,-Q both in E(K), of which exactly one is a
            # solution.  So we just check the existence of roots:
            return len(g.roots())>0

        # Now 2*P==0

        if m%2==1:
            return True # P itself is a solution when m is odd

        # Now m is even and 2*P=0.  Roots of g in K may or may not
        # lift to solutions in E(K), so we fall back to the default.
        # Note that division polynomials are cached so this is not
        # inefficient:

        return len(self.division_points(m)) > 0

    def division_points(self, m, poly_only=False):
        r"""
        Return a list of all points `Q` such that `mQ=P` where `P` = self.

        Only points on the elliptic curve containing self and defined
        over the base field are included.

        INPUT:

        - ``m`` -- a positive integer
        - ``poly_only`` -- bool (default: False); if True return polynomial whose roots give all possible `x`-coordinates of `m`-th roots of self.

        OUTPUT:

        (list) -- a (possibly empty) list of solutions `Q` to `mQ=P`,  where `P` = self.

        EXAMPLES:

        We find the five 5-torsion points on an elliptic curve::

            sage: E = EllipticCurve('11a'); E
            Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
            sage: P = E(0); P
            (0 : 1 : 0)
            sage: P.division_points(5)
            [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)]

        Note above that 0 is included since [5]*0 = 0.

        We create a curve of rank 1 with no torsion and do a consistency check::

            sage: E = EllipticCurve('11a').quadratic_twist(-7)
            sage: Q = E([44,-270])
            sage: (4*Q).division_points(4)
            [(44 : -270 : 1)]

        We create a curve over a non-prime finite field with group of order `18`::

            sage: k.<a> = GF(25)
            sage: E = EllipticCurve(k, [1,2+a,3,4*a,2])
            sage: P = E([3,3*a+4])
            sage: factor(E.order())
            2 * 3^2
            sage: P.order()
            9

        We find the `1`-division points as a consistency check -- there
        is just one, of course::

            sage: P.division_points(1)
            [(3 : 3*a + 4 : 1)]

        The point `P` has order coprime to 2 but divisible by 3, so::

            sage: P.division_points(2)
            [(2*a + 1 : 3*a + 4 : 1), (3*a + 1 : a : 1)]

        We check that each of the 2-division points works as claimed::

            sage: [2*Q for Q in P.division_points(2)]
            [(3 : 3*a + 4 : 1), (3 : 3*a + 4 : 1)]

        Some other checks::

            sage: P.division_points(3)
            []
            sage: P.division_points(4)
            [(0 : 3*a + 2 : 1), (1 : 0 : 1)]
            sage: P.division_points(5)
            [(1 : 1 : 1)]

        An example over a number field (see trac #3383)::

            sage: E = EllipticCurve('19a1')
            sage: K.<t> = NumberField(x^9-3*x^8-4*x^7+16*x^6-3*x^5-21*x^4+5*x^3+7*x^2-7*x+1)
            sage: EK = E.base_extend(K)
            sage: E(0).division_points(3)
            [(0 : 1 : 0), (5 : -10 : 1), (5 : 9 : 1)]
            sage: EK(0).division_points(3)
            [(0 : 1 : 0), (5 : 9 : 1), (5 : -10 : 1)]
            sage: E(0).division_points(9) 
            [(0 : 1 : 0), (5 : -10 : 1), (5 : 9 : 1)]
            sage: EK(0).division_points(9)
            [(0 : 1 : 0), (5 : 9 : 1), (5 : -10 : 1), (-150/121*t^8 + 414/121*t^7 + 1481/242*t^6 - 2382/121*t^5 - 103/242*t^4 + 629/22*t^3 - 367/242*t^2 - 1307/121*t + 625/121 : 35/484*t^8 - 133/242*t^7 + 445/242*t^6 - 799/242*t^5 + 373/484*t^4 + 113/22*t^3 - 2355/484*t^2 - 753/242*t + 1165/484 : 1), (-150/121*t^8 + 414/121*t^7 + 1481/242*t^6 - 2382/121*t^5 - 103/242*t^4 + 629/22*t^3 - 367/242*t^2 - 1307/121*t + 625/121 : -35/484*t^8 + 133/242*t^7 - 445/242*t^6 + 799/242*t^5 - 373/484*t^4 - 113/22*t^3 + 2355/484*t^2 + 753/242*t - 1649/484 : 1), (-1383/484*t^8 + 970/121*t^7 + 3159/242*t^6 - 5211/121*t^5 + 37/484*t^4 + 654/11*t^3 - 909/484*t^2 - 4831/242*t + 6791/484 : 927/121*t^8 - 5209/242*t^7 - 8187/242*t^6 + 27975/242*t^5 - 1147/242*t^4 - 1729/11*t^3 + 1566/121*t^2 + 12873/242*t - 10871/242 : 1), (-1383/484*t^8 + 970/121*t^7 + 3159/242*t^6 - 5211/121*t^5 + 37/484*t^4 + 654/11*t^3 - 909/484*t^2 - 4831/242*t + 6791/484 : -927/121*t^8 + 5209/242*t^7 + 8187/242*t^6 - 27975/242*t^5 + 1147/242*t^4 + 1729/11*t^3 - 1566/121*t^2 - 12873/242*t + 10629/242 : 1), (-4793/484*t^8 + 6791/242*t^7 + 10727/242*t^6 - 18301/121*t^5 + 2347/484*t^4 + 2293/11*t^3 - 7311/484*t^2 - 17239/242*t + 26767/484 : 30847/484*t^8 - 21789/121*t^7 - 34605/121*t^6 + 117164/121*t^5 - 10633/484*t^4 - 29437/22*t^3 + 39725/484*t^2 + 55428/121*t - 176909/484 : 1), (-4793/484*t^8 + 6791/242*t^7 + 10727/242*t^6 - 18301/121*t^5 + 2347/484*t^4 + 2293/11*t^3 - 7311/484*t^2 - 17239/242*t + 26767/484 : -30847/484*t^8 + 21789/121*t^7 + 34605/121*t^6 - 117164/121*t^5 + 10633/484*t^4 + 29437/22*t^3 - 39725/484*t^2 - 55428/121*t + 176425/484 : 1)]

        """
        # Coerce the input m to an integer
        m = rings.Integer(m)
        # Check for trivial cases of m = 1, -1 and 0.
        if m == 1 or m == -1:
            return [self]
        if m == 0:
            if self == 0: # then every point Q is a solution, but...
                return [self]
            else:
                return []

        # ans will contain the list of division points. 
        ans = []

        # We compute a polynomial g whose roots give all possible x
        # coordinates of the m-division points.  The number of
        # solutions (over the algebraic closure) is m^2, assuming that
        # the characteristic does not divide m.

        E = self.curve()
        P = self
        nP = -P
        P_is_2_torsion = (P==nP)

        # If self is the 0, then self is a solution, and the correct
        # poly is the m'th division polynomial
        if P == 0:
            ans.append(P)
            g = E.division_polynomial(m)
        else:
            # The poly g here is 0 at x(Q) iff x(m*Q) = x(P).
            g = E._multiple_x_numerator(m) - P[0]*E._multiple_x_denominator(m)

            # When 2*P=0, then -Q is a solution iff Q is.  For even m,
            # no 2-torsion point is a solution, so that g is the
            # square of a polynomial g1 of degree m^2/2, and each root
            # of g1 leads to a pair of solutions Q, -Q to m*Q=P.  For
            # odd m, P itself is the only 2-torsion solution, so g has
            # the form (x-x(P))*g1(x)^2 where g1 has degree (m^2-1)/2
            # and each root of g1 leads to a pair Q, -Q.

            if P_is_2_torsion:
                if m%2==0:
                    # This computes g.sqrt() which is not implemented
                    g = g.gcd(g.derivative())*g.leading_coefficient().sqrt()

            # When 2*P!=0, then for each solution Q to m*Q=P, -Q is
            # not a solution (and points of order 2 are not
            # solutions).  Hence the roots of g are distinct and each
            # gives rise to precisely one solution Q.
            
                else:
                    g0 = g.variables()[0] - P[0]
                    g = g // g0
                    g = g.gcd(g.derivative())*g.leading_coefficient().sqrt()
                    g = g0*g
                    
        if poly_only:
            return g

        for x in g.roots(multiplicities=False):
            if E.is_x_coord(x):
                # Make a point on the curve with this x coordinate.
                Q = E.lift_x(x)
                nQ = -Q
                mQ = m*Q
                # if P==-P then Q works iff -Q works, so we include
                # both unless they are equal:
                if P_is_2_torsion:
                    if mQ == P:
                        ans.append(Q)
                        if nQ != Q:
                            ans.append(nQ)
                else:
                    # P is not 2-torsion so at most one of Q, -Q works
                    # and we must try both:
                    if mQ == P:
                        ans.append(Q)
                    elif mQ == nP:
                        ans.append(nQ)
                        
        # Finally, sort and return
        ans.sort()
        return ans

    def _divide_out(self,p):
        r"""
        Return `(Q,k)` where `p^kQ` == self and `Q` cannot be divided by `p`.

        ..WARNING: 

        It is up to the caller to make sure that this does not loop
        endlessly.  It is used in
        ``EllipticCurve_generic._p_primary_torsion_basis()``, when
        self will always have (finite) order which is a power of `p`,
        so that the order of `Q` increases by a factor of `p` at each
        stage.

        Since it will clearly be in danger of looping when
        self.is_zero(), this case is caught, but otherwise caveat
        user.

        EXAMPLES::

            sage: E = EllipticCurve('37a1')
            sage: P = E([0, 0])
            sage: R = 12*P
            sage: R._divide_out(2)
            ((-1 : -1 : 1), 2)
            sage: R._divide_out(3)
            ((2 : -3 : 1), 1)
            sage: R._divide_out(5)
            ((1357/841 : 28888/24389 : 1), 0)
            sage: R._divide_out(12)
            Traceback (most recent call last):
            ...
            ValueError: p (=12) should be prime.
        """
        p = rings.Integer(p)
        if not p.is_prime():
            raise ValueError, "p (=%s) should be prime."%p

        if self.is_zero():
            raise ValueError, "self must not be 0."

        k=0; Q=self
        pts = Q.division_points(p)
        while len(pts) > 0:
            Q = pts[0]
            k += 1
            pts = Q.division_points(p)
        return (Q,k)


    ##############################  end  ################################

    def _line_(self,R,Q):
        r""" 
        Computes the value at `Q` of a straight line through points self and `R`.

        INPUT:

        - ``R, Q`` -- points on self.curve() with ``Q`` nonzero.

        OUTPUT:

        An element of the base field self.curve().base_field().
        A ValueError is raised if ``Q`` is zero.

        EXAMPLES::

            sage: F.<a>=GF(2^5)
            sage: E=EllipticCurve(F,[0,0,1,1,1])
            sage: P = E(a^4 + 1, a^3)
            sage: Q = E(a^4, a^4 + a^3)
            sage: O = E(0)
            sage: P._line_(P,-2*P) == 0
            True
            sage: P._line_(Q,-(P+Q)) == 0
            True
            sage: O._line_(O,Q) == F(1)
            True
            sage: P._line_(O,Q) == a^4 - a^4 + 1 
            True
            sage: P._line_(13*P,Q) == a^4 
            True
            sage: P._line_(P,Q) == a^4 + a^3 + a^2 + 1
            True

        See trac #7116::
            
            sage: P._line_ (Q,O)
            Traceback (most recent call last):
            ...
            ValueError: Q must be nonzero.

        ..NOTES:  

            This function is used in _miller_ algorithm.

        AUTHOR:

        - David Hansen (2009-01-25)
        """
        if Q.is_zero():
            raise ValueError, "Q must be nonzero."

        if self.is_zero() or R.is_zero():
            if self == R:
                return self.curve().base_field().one_element()
            if self.is_zero():
                return Q[0] - R[0]
            if R.is_zero():
                return Q[0] - self[0]
        elif self != R:
            if self[0] == R[0]:
                return Q[0] - self[0]
            else:
                l = (R[1] - self[1])/(R[0] - self[0])
                return Q[1] - self[1] - l * (Q[0] - self[0])
        else:
            a1, a2, a3, a4, a6 = self.curve().a_invariants()
            numerator = (3*self[0]**2 + 2*a2*self[0] + a4 - a1*self[1])
            denominator = (2*self[1] + a1*self[0] + a3)
            if denominator == 0:
                return Q[0] - self[0]
            else:
                l = numerator/denominator
                return Q[1] - self[1] - l * (Q[0] - self[0])

    def _miller_(self,Q,n):
        r"""
        Compute the value at `Q` of the rational function `f_{n,P}`, where the divisor of `f_{n,P}` is `n[P]-n[O]`.

        INPUT:

        - ``Q`` -- a nonzero point on self.curve().

        - ``n`` -- an integer such that `n*P = n*Q = (0:1:0)` where `P`=self.

        OUTPUT:

        An element in the base field self.curve().base_field()
        A ValueError is raised if ``Q`` is zero.        

        EXAMPLES::

            sage: F.<a>=GF(2^5)
            sage: E=EllipticCurve(F,[0,0,1,1,1])
            sage: P = E(a^4 + 1, a^3)
            sage: Fx.<b>=GF(2^(4*5))
            sage: Ex=EllipticCurve(Fx,[0,0,1,1,1])
            sage: phi=Hom(F,Fx)(F.gen().minpoly().roots(Fx)[0][0])
            sage: Px=Ex(phi(P.xy()[0]),phi(P.xy()[1]))
            sage: Qx = Ex(b^19 + b^18 + b^16 + b^12 + b^10 + b^9 + b^8 + b^5 + b^3 + 1, b^18 + b^13 + b^10 + b^8 + b^5 + b^4 + b^3 + b) 
            sage: Px._miller_(Qx,41) == b^17 + b^13 + b^12 + b^9 + b^8 + b^6 + b^4 + 1
            True
            sage: Qx._miller_(Px,41) == b^13 + b^10 + b^8 + b^7 + b^6 + b^5
            True
            sage: P._miller_(E(0),41)
            Traceback (most recent call last):
            ...
            ValueError: Q must be nonzero.

        An example of even order::

            sage: F.<a> = GF(19^4)
            sage: E = EllipticCurve(F,[-1,0])
            sage: P = E(15*a^3 + 17*a^2 + 14*a + 13,16*a^3 + 7*a^2 + a + 18)
            sage: Q = E(10*a^3 + 16*a^2 + 4*a + 2, 6*a^3 + 4*a^2 + 3*a + 2)
            sage: x=P.weil_pairing(Q,360)
            sage: x^360 == F(1)
            True

        You can use the _miller_ function on linearly dependent points, but with the risk of a dividing with zero::

            sage: Px._miller_(2*Px,41)
            Traceback (most recent call last):
            ...
            ZeroDivisionError: division by zero in finite field.

        ALGORITHM: 

            Double-and-add.

        REFERENCES:

        - [Mil04] Victor S. Miller, "The Weil pairing, and its efficient calculation", J. Cryptol., 17(4):235-261, 2004

        AUTHOR:

        - David Hansen (2009-01-25)

        """
        if Q.is_zero():
            raise ValueError, "Q must be nonzero."

        t = self.curve().base_field().one_element()
        V = self
        S = 2*V
        nbin = n.bits()
        i = n.nbits() - 2
        while i > -1:
            S = 2*V
            t = (t**2)*(V._line_(V,Q)/S._line_(-S,Q)) 
            V = S
            if nbin[i] == 1:
                S = V+self
                t=t*(V._line_(self,Q)/S._line_(-S,Q))
                V = S
            i=i-1
        return t

    def weil_pairing(self, Q, n):
        r"""
        Compute the Weil pairing of self and `Q` using Miller's algorithm.

        INPUT:

        - ``Q`` -- a point on self.curve().

        - ``n`` -- an integer `n` such that `nP = nQ = (0:1:0)` where `P` = self.

        OUTPUT:

        An `n`'th root of unity in the base field self.curve().base_field()

        EXAMPLES::

            sage: F.<a>=GF(2^5)
            sage: E=EllipticCurve(F,[0,0,1,1,1])
            sage: P = E(a^4 + 1, a^3)
            sage: Fx.<b>=GF(2^(4*5))
            sage: Ex=EllipticCurve(Fx,[0,0,1,1,1])
            sage: phi=Hom(F,Fx)(F.gen().minpoly().roots(Fx)[0][0])
            sage: Px=Ex(phi(P.xy()[0]),phi(P.xy()[1]))
            sage: O = Ex(0)
            sage: Qx = Ex(b^19 + b^18 + b^16 + b^12 + b^10 + b^9 + b^8 + b^5 + b^3 + 1, b^18 + b^13 + b^10 + b^8 + b^5 + b^4 + b^3 + b) 
            sage: Px.weil_pairing(Qx,41) == b^19 + b^15 + b^9 + b^8 + b^6 + b^4 + b^3 + b^2 + 1
            True
            sage: Px.weil_pairing(17*Px,41) == Fx(1)
            True
            sage: Px.weil_pairing(O,41) == Fx(1)
            True

        An error is raised if either point is not n-torsion::

            sage: Px.weil_pairing(O,40)
            Traceback (most recent call last):
            ...
            ValueError: points must both be n-torsion

        A larger example (see trac \#4964)::

            sage: P,Q = EllipticCurve(GF(19^4,'a'),[-1,0]).gens()
            sage: P.order(), Q.order()
            (360, 360)
            sage: z = P.weil_pairing(Q,360)
            sage: z.multiplicative_order()
            360

        An example over a number field::

            sage: P,Q = EllipticCurve('11a1').change_ring(CyclotomicField(5)).torsion_subgroup().gens()  # long time (10s)
            sage: P,Q = (P.element(), Q.element())  # long time
            sage: (P.order(),Q.order())  # long time
            (5, 5)
            sage: P.weil_pairing(Q,5)  # long time
            zeta5^2
            sage: Q.weil_pairing(P,5)  # long time
            zeta5^3        

        ALGORITHM: 

        Implemented using Proposition 8 in [Mil04].  The value 1 is
        returned for linearly dependent input points.  This condition
        is caught via a DivisionByZeroError, since the use of a
        discrete logarithm test for linear dependence, is much to slow
        for large `n`.

        REFERENCES: 

        [Mil04] Victor S. Miller, "The Weil pairing, and its efficient
        calculation", J. Cryptol., 17(4):235-261, 2004

        AUTHOR:

        - David Hansen (2009-01-25)
        """
        P = self
        E = P.curve()
        
        if not Q.curve() is E:
            raise ValueError, "points must both be on the same curve"
            
        # Test if P, Q are both in E[n]
        if not ((n*P).is_zero() and (n*Q).is_zero()):
            raise ValueError, "points must both be n-torsion"

        one = E.base_field().one_element()

        # Case where P = Q
        if P == Q:
            return one

        # Case where P = O or Q = O
        if P.is_zero() or Q.is_zero():
            return one

        # The non-trivial case P != Q

        # Note (2010-12-29): The following code block should not be
        # used.  It attempts to reduce the pairing calculation to order
        # d = gcd(|P|,|Q|) instead of order n, but this approach is
        # counterproductive, since calculating |P| and |Q| is much
        # slower than calculating the pairing [BGN05].
        # 
        # [BGN05] D. Boneh, E. Goh, and K. Nissim, "Evaluating 2-DNF
        # Formulas on Ciphertexts", TCC 2005, LNCS 3378, pp. 325-341.
        #
        # Reduction to order d = gcd(|P|,|Q|); value is a d'th root of unity
        # try:
        #     nP = P.order()
        # except AttributeError:
        #     nP = generic.order_from_multiple(P,n,operation='+')
        # try:
        #     nQ = Q.order()
        # except AttributeError:
        #     nQ = generic.order_from_multiple(Q,n,operation='+')
        # d = arith.gcd(nP,nQ)
        # if d==1:
        #     return one
        #
        # P = (nP//d)*P # order d
        # Q = (nQ//d)*Q # order d
        # n = d

        try: 
            return ((-1)**n.test_bit(0))*(P._miller_(Q,n)/Q._miller_(P,n))
        except ZeroDivisionError, detail:
            return one

class EllipticCurvePoint_number_field(EllipticCurvePoint_field):
    """
    A point on an elliptic curve over a number field.

    Most of the functionality is derived from the parent class
    ``EllipticCurvePoint_field``.  In addition we have support for the
    order of a point, and heights (currently only implemented over
    `\QQ`).

    EXAMPLES::

        sage: E = EllipticCurve('37a')
        sage: E([0,0])
        (0 : 0 : 1)
        sage: E(0,0)               # brackets are optional
        (0 : 0 : 1)
        sage: E([GF(5)(0), 0])     # entries are coerced
        (0 : 0 : 1)      

        sage: E(0.000, 0)
        (0 : 0 : 1)
         
        sage: E(1,0,0)  
        Traceback (most recent call last):
        ...
        TypeError: Coordinates [1, 0, 0] do not define a point on
        Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field

    ::

        sage: E = EllipticCurve([0,0,1,-1,0])
        sage: S = E(QQ); S
        Abelian group of points on Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field

    TESTS::

        sage: loads(S.dumps()) == S
        True
        sage: P = E(0,0); P
        (0 : 0 : 1)
        sage: loads(P.dumps()) == P
        True
        sage: T = 100*P
        sage: loads(T.dumps()) == T
        True

    Test pickling an elliptic curve that has known points on it::

        sage: e = EllipticCurve([0, 0, 1, -1, 0]); g = e.gens(); loads(dumps(e)) == e
        True
    """

    def order(self):
        r"""
        Return the order of this point on the elliptic curve.

        If the point has infinite order, returns +Infinity.  For
        curves defined over `\QQ`, we call PARI; over other
        number fields we implement the function here.

        .. note::

           :meth:`additive_order` is a synonym for :meth:`order`

        EXAMPLES::

            sage: E = EllipticCurve([0,0,1,-1,0])
            sage: P = E([0,0]); P
            (0 : 0 : 1)
            sage: P.order()
            +Infinity

        ::

            sage: E = EllipticCurve([0,1])
            sage: P = E([-1,0])
            sage: P.order()
            2
            sage: P.additive_order()
            2

        """
        try:
            return self._order
        except AttributeError:
            pass

        if self.is_zero():
            self._order = rings.Integer(1)
            return self._order

        E = self.curve()

        # Special code for curves over Q, calling PARI
        try:
            n = int(E.pari_curve().ellorder(self))
            if n == 0: n = oo
            self._order = n
            return n
        except PariError:
            pass

        # Get the torsion order if known, else a bound on (multiple
        # of) the order.  We do not compute the torsion if it is not
        # already known, since computing the bound is faster (and is
        # also cached).

        try:
            N = E._torsion_order
        except AttributeError:
            N = E._torsion_bound()

        # Now self is a torsion point iff it is killed by N:
        if not (N*self).is_zero():
            self._order = oo
            return self._order
            
        # Finally we find the exact order using the generic code:
        self._order = generic.order_from_multiple(self,N,operation='+')
        return self._order

    additive_order = order
    
    def has_finite_order(self):
        """
        Return True iff this point has finite order on the elliptic curve.

        EXAMPLES::

            sage: E = EllipticCurve([0,0,1,-1,0])
            sage: P = E([0,0]); P
            (0 : 0 : 1)
            sage: P.has_finite_order()
            False

        ::

            sage: E = EllipticCurve([0,1])
            sage: P = E([-1,0])
            sa…