class LTorsionGroup( metaclass=template("_elliptic_curve") ):
    """
    An l-torsion group of an elliptic curve for odd l; use the @c elements()
    method to have intuitive syntax when working with l-torsion points.
    
    This is a template class that must be instantiated with the elliptic curve.
    Use it, for example, as follows:
    @code
    # Create the underlying elliptic curve
    E = elliptic_curves.naive.EllipticCurve( FiniteField(7), 3, 4 )
    # Instantiate the template; El is a class (here: l-torsion groups of E)
    El = LTorsionGroup( E )
    
    # Construct the group of 3-torsion points and do something with its points
    for P in El[3].elements():
        do_something(P)
    @endcode
    
    The l-Torsion points of an elliptic curve are the points of order l.  Their
    coordinates come from the algebraic closure of the elliptic curve's field,
    which makes them badly suited for computer representation. (Extending the
    coordinate domain to the algebraic closure makes sure that all points of
    order l are available; there are always @f$ l^2 @f$ points of order l if
    l is prime to the field characteristic.) 
    
    @note  The l-torsion points are represented implicitly through polynomial
           arithmetic.  Therefore, the list of points in returned by
           @c elements() contains only a single entry: the point @f$ (x, y) \in
           E\bigl(\mathbb{F}_{p}(E) / \psi_l\bigr) @f$, where
           @f$ \mathbb{F}_{p}(E)/\psi_l @f$ denotes the field of rational
           functions over the elliptic curve @f$ E @f$ with numerator and
           denominator polynomials taken modulo the l-th division polynomial
           @f$ \psi_l @f$.
           
    @note  The class supports only odd torsions greater one because
           multivariate polynomial arithmetic is unavailable
           
    @see   Silverman, Joseph H., "The Arithmetic of Elliptic Curves",
           second edition, Springer, 2009, p. 373
    """

    #- Instance Methods ------------------------------------------------------- 
    
    def __init__(self, torsion):
        """
        Construct a new l-torsion group for the given @p torsion.
        
        @param torsion The @p torsion is the order of the points in the group;
                       the group represents all points of order @p torsion.  It
                       must be an odd integer greater than 1 and prime to the
                       field characteristic; the limitation comes from the
                       implicit representation, see above. 
        """
        torsion = int( torsion )
        if torsion < 3 or torsion % 2 == 0 or self.curve().field()(torsion) == 0:
            raise ValueError( "only odd torsions greater than 1 are supported" )
        
        self.__torsion = torsion
        self.__point = None
    
    
    def elements(self):
        """
        Return a list containing the one point that implicitly represents
        the whole group.
        
        Use it, for example, to elegantly perform an operation on all l-torsion
        points:
        @code
        for P in torsion_group.elements():
            do_something(P)
        @endcode

        The one point in the list is @f$ (x, y) \in
        E\bigl(\mathbb{F}_{p}(E) / \psi_l\bigr) @f$, where
        @f$ \mathbb{F}_{p}(E)/\psi_l @f$ denotes the field of rational functions
        over the elliptic curve @f$ E @f$ with numerator and denominator
        polynomials taken modulo the l-th division polynomial @f$ \psi_l @f$.
        """
        if not self.__point:
            self.__init_point()
        return [ self.__point ]
    
    
    def torsion(self):
        """
        Return the group's torsion, that is, l for the l-torsion group.
        """
        return self.__torsion


    def __init_point(self):
        """
        Create the point that implicitly represents the whole l-torsion group
        and store it in the cache.
        
        Calling the method has no effect if the point already exists.
        """
        if self.__point:
            return
        
        # R = F[x] / (y**2 - x**3 - A*x - B)
        R = self._division_polynomial_list().curve_polynomials()
        # The n-th division polynomial
        psi = self._division_polynomial_list()[ self.__torsion ]

        # T = ( F[x] / (y**2 - x**3 - A*x - B) ) / psi(l)
        S = QuotientRing( R, psi )
        T = FractionField( S )
        
        A, B = R.curve().parameters()
        
        # Polynomials x and y on the curve interpreted
        # as elements of the field of fractions
        x = T( R( (0, 1), 0 ) )
        y = T( R( 0     , 1 ) )
        
        self.__point = EllipticCurve( T, A, B )( x, y )


    #- Class Methods----------------------------------------------------------- 
    
    @classmethod
    def curve(cls):
        """
        Return the elliptic curve that was used to initialize the template.
        """
        return cls._elliptic_curve
    
    @classmethod
    def _division_polynomial_list(cls):
        """
        Return the list of division polynomials that supplies the modulus for
        the implicit polynomial representation of the l-torsion points.
        """
        try:
            return cls.__division_polynomial_list
        except AttributeError:
            # R = F[x] / (y**2 - x**3 - A*x - B)
            # where F, A, B are the elliptic curve's field and its parameters.
            R = CurvePolynomials( cls._elliptic_curve )
            cls.__division_polynomial_list = DivisionPolynomialsList( R )

            return cls.__division_polynomial_list def inverse_primorial(n, shunned = 0):
    """
    Return the smallest prime @c p such that the product of all primes
    up to @c p (inclusive) is at least @p n. For example,
    @c inverse_primorial(30) is 5, and @c inverse_primorial(31) is 7. 
    
    @note      This function uses primes_range() to obtain a list of primes.
               See the notes there for use case limitations.

    @param     n   The number object that the product must exceed in size.
                   The object must have an integer interpretation.
    
    @return    The prime @c p such that the product of all primes up to @c p
               (inclusive) is at least @p n; if @p n is too small (less than 2)
               the result is 2.
    """
    product = 1
    # A much smaller upper bound should suffice; however, we play it safe.
    for prime in primes_range(2, int(n)+1):
        if prime != shunned:
            product *= prime
            if product >= n:
                return prime
    # Return the smallest prime if n is too small
    return 2 def primes_range(lower_bound, upper_bound):
    """
    Return a list of all primes in the range [@p lower_bound, @p upper_bound).
    
    @note      The function returns an actual list, not a proper range object.
               Furthermore, it computes the complete list of primes up to
               @p upper_bound on each call. These characteristics make it badly
               suited for anything involving large ranges or large primes.

    @param     lower_bound     The minimum size for returned primes; typically
                               this is an integer, but any number object with
                               an integer interpretation works. 
    @param     upper_bound     The strict upper bound for returned primes: all
                               returned primes are strictly less than this
                               number.
    """
    # The sieve of Eratosthenes
    primes = set( range(2, upper_bound) )
    sieve_bound = ceil( sqrt(upper_bound) )
    for i in range(2, sieve_bound):
        if i in primes:
            multiples = [ i*j for j in range(2, ceil(upper_bound / i)) ]
            primes -= set(multiples)
    return sorted(list( primes - set(range(2, lower_bound)) )) def solve_congruence_equations( congruences ):
    """
    Return a quotient class that simultaneously solves the list of
    @p congruences; this is the Chinese Remainder Theorem.
    
    All representatives of the returned congruence have the remainders of
    @p congruences when taken modulo the respective moduli. The result is a
    congruence modulo the product of all moduli. Thus the function returns a
    number @f$ z \mod \prod_{i} m_i @f$ such that
    @f{align*}{
     z &\equiv a_1 \mod m_1 \\
     \vdots \\
     z &\equiv a_k \mod m_k
    @f}
    
    @note The moduli @f$ m_i @f$ of all congruences must be relatively prime.
    
    @exception ValueError      if @p congruences is empty.

    @param     congruences     An iterable of objects of QuotientClass over the
                               Integers. Every pair of two different moduli
                               must have a greatest common divisor (gcd()) of 1.
    
    @return    An instance of rings.quotients.naive.QuotientRing over the
               rings.integers.naive.Integers solving the @p congruences.
    
    @see       Robinson, D. J. S., "An Introduction to Abstract Algebra", p. 27
    """
    # The Chinese remainder theorem
    if not congruences:
        raise ValueError( "cannot solve empty equation system" )
        
    if __debug__:
        # This test is expensive; remove it in optimized execution
        pairs = [ (c1, c2) for c1 in congruences for c2 in congruences if c1 != c2 ]
        pairwise_gcds = [ gcd( c1.modulus(), c2.modulus() ) for c1, c2 in pairs ]
    
        assert set( pairwise_gcds ) == set([ 1 ]), \
            "the Chinese Remainder Theorem requires relatively prime moduli"
    
    common_modulus = reduce( mul, [ c.modulus() for c in congruences ] )
    common_representative = 0
    for c in congruences:
        neutralizer = common_modulus // c.modulus()
        common_representative += c.remainder() * neutralizer  \
                                    * inverse_modulo( neutralizer, c.modulus() )
    
    quotient_ring = QuotientRing( Integers, common_modulus )
    return quotient_ring( common_representative ) def representative_in_range( quotient_class, valid_range ):
    """
    Return the unique representative of @p quotient_class in @p valid_range.
    For example, the representative of @c (3 mod 7) in @c range(13, 19) is 17.
    
    @param     quotient_class  A QuotientClass over the Integers.
    @param     valid_range     The range() object describing the interval in
                               which to the returned representative must lie. 
    
    @return    An integer @c x in @p valid_range with @c x % m == r, where
               @c m is the modulus, and @c r the remainder of @p quotient_class.  
    
    @exception ValueError  either if there is no representative in
                           @p valid_range, or if there are multiple.
    """
    if len(valid_range) > quotient_class.modulus():
        raise ValueError("solution not unique")
    
    # range[0] is the first element, that is, the lower bound
    q, r = divmod( valid_range[0], quotient_class.modulus() )
    shifted_range = range(r, r + len(valid_range))
    
    if quotient_class.remainder() in shifted_range:
        return q * quotient_class.modulus()  + quotient_class.remainder()
    elif quotient_class.remainder() + quotient_class.modulus() in shifted_range:
        return (q+1) * quotient_class.modulus()  + quotient_class.remainder()
    else:
        # FIXME: Really use exceptions?
        raise ValueError("no solution") class Integers(int):
    """
    The ring of integers.  This wrapper class extends Python's built-in integer
    type @c int with the interface provided by the other ring classes.
    
    @see   rings.polynomials.naive.Polynomials,
           rings.quotients.naive.QuotientRing
    """
    
    @staticmethod
    def zero():
        """
        Return @c 0, the neutral element of integer addition.
        """
        return 0

    @staticmethod
    def one():
        """
        Return @c 1, the neutral element of integer multiplication.
        """
        return 1 class Polynomials( CommutativeRing, metaclass=template( "_coefficient_field" ) ):
    """
    A ring of polynomials in one indeterminate with coefficients from the
    provided field; the polynomials support infix notation for ring operations,
    as well as mixed argument types.
    
    This is a template class that must be instantiated with the coefficient
    field.  Use it, for example, as follows:
    @code
    # Instantiate the template; GF7x is a class (here: polynomials over GF7)
    GF7 = fields.finite.naive.FiniteField( 7 )
    GF7x = Polynomials( GF7 )
    p = GF7x(1, 1)    # Create a polynomial: x + 1
    q = GF7x(-1, 1)   # Another polynomial: x - 1
    p * q == GF7x(-1, 0, 1)    # This is true: (x + 1)(x - 1) = x^2 - 1
    type(p) is GF7x   # This is also true

    # Use the polynomial 'x' for a more natural syntax
    x = GF7(0, 1)
    x**3 - 2*x + 1 == GF7x(1, -2, 0, 1)
    @endcode
    
    The polynomial ring operations are defined in terms of the coefficient
    field operations.  For everything to work as expected, the template
    argument must be a class whose instances implement the fields.Field
    interface.  Typically, this will be some specialization of
    fields.finite.naive.FiniteField.
    
    Rings, however, also work as a source of coefficients.  Only division
    requires that the leading coefficient is a unit. 

    @note  The implementation emphasizes simplicity over speed; it omits
           possible optimizations.
    
    @note  The class uses the operand_casting() decorator: @c other operands in
           binary operations will first be treated as Polynomial elements.
           If that fails, the operation will be repeated after @p other was fed
           to the constructor __init__().  If that fails, too, then the
           operation returns @c NotImplemented (so that @p other.__rop__()
           might be invoked).
    
    @see   For example, Robinson, Derek J. S.,
           "An Introduction to Abstract Algebra", p. 100.
    """    

    #- Instance Methods ----------------------------------------------------------- 
    
    def __init__(self, element_description, *further_coefficients):
        """
        Create a new polynomial from the given @p element_description.
        
        @param element_description
            Valid input values are:
            - A Polynomial over the same field: then the polynomial is copied
              and @p further_coefficients is ignored.
            - An iterable object (with an @c __iter__() method): then the
              polynomial coefficients are read from the iterable and, again,
              @p further_coefficients is ignored.  The list must be in
              ascending order: first the constant, then the linear, quadratic,
              and cubic coefficients; and so on.
            - Any number of arguments that can be interpreted as coefficient
              field elements; they will be combined
              into the list of coefficients
        
        @param further_coefficients    The list of variable positional
                                       arguments.
        
        For example, there are several ways to construct the polynomial
        @f$ x^3 + 4x^2 - 3x + 2@f$:
        @code
        # Suppose R is a Polynomials template specialization
        # List of coefficients
        p2 = R( [2, -3, 4, 1] )
        # Variable argument count
        p3 = R( 2, -3, 4, 1 )
        # Copy
        p1 = R( p2 )
        @endcode
        
        @note  Leading zeros will be ignored.
        """
        # List of coefficients is in ascending order without leading zeros.
        # Example: x^2 + 2x + 5 = [5, 2, 1]
        if isinstance(element_description, self.__class__):
            # A reference suffices because objects are immutable
            self.__coefficients = element_description.__coefficients
        
        else:
            if type(element_description) in [ list, tuple ]:
                coefficients = element_description
            else:
                coefficients = [ element_description ] + list(further_coefficients)  

            F = self._coefficient_field
            self.__coefficients = [ F(c) for c in coefficients ]

        self.__remove_leading_zeros()

    
    def coefficients(self):
        """
        Return a copy of the list of coefficients.  The list is in ascending
        order: it starts with the constant coefficient, then the linear,
        quadratic, and cubic coefficients; and so on.  It ends with the leading
        coefficient.
        """
        return self.__coefficients[:]


    def leading_coefficient(self):
        """
        Return the leading coefficient, or zero if @p self is
        the zero polynomial.
        """
        if self.__coefficients:
            return self.__coefficients[-1]
        else:
            return self._coefficient_field.zero()


    def degree(self):
        """
        Return the degree of the polynomial.
        
        @note  The zero polynomial has degree @f$ -\infty @f$; however, this
               implementation returns the number @f$ -2^{30} @f$ to avoid
               special constructions that side-step Python's @c int objects.
               For practical purposes, this should suffice.
        """
        if not self.__coefficients:
            # FIXME: The degree of the zero polynomial is minus infinity.
            #        This will, however, do for now.
            return -( 2**30 )
        return len( self.__coefficients ) - 1

    
    def __bool__(self):
        """
        Test whether the polynomial is non-zero: return @c True if, and only
        if, at least one coefficient is non-zero. Return @c False if all
        coefficients are zero.
        
        Implicit conversions to boolean (truth) values use this method, for
        example when @c x is an element of Polynomials:
        @code
        if x:
            do_something()
        @endcode
        """
        return bool( self.__coefficients )


    def __eq__(self, other):
        """
        Test whether another polynomial @p other is equal to @p self; return
        @c True if that is the case.  The infix operator @c == calls
        this method.
        
        Two polynomials are equal if, and only if, all their coefficients are
        equal. 
        """
        if self.degree() != other.degree():
            return False
        
        # zip() is OK because we have identical length
        for x, y in zip(self.__coefficients, other.__coefficients):
            if x != y:
                return False

        return True


    def __add__(self, other):
        """
        Return the sum of @p self and @p other. The infix operator @c + calls
        this method.
        
        Polynomials add coefficient-wise without carrying: the sum of two
        polynomials @f$ \sum_{k} a_{k}x^{k} @f$ and
        @f$ \sum_{k} b_{k}x^{k} @f$ is the polynomial
        @f$ \sum_{k} (a_{k} + b_{k})x^{k} @f$.
        """
        zero = self._coefficient_field.zero()
        coefficient_pairs = self.__pad_and_zip(
                                    self.__coefficients,
                                    other.__coefficients,
                                    zero
                                )
        coefficient_sums = [ x + y for x, y in coefficient_pairs ]
        return self.__class__( coefficient_sums )
    
    
    def __neg__(self):
        """
        Return the additive inverse (the negative) of @p self.
        
        A polynomial is negated by negating its coefficients: the additive
        inverse of @f$ \sum_{k} a_{k}x^{k} @f$ is @f$ \sum_{k} -a_{k}x^{k} @f$.
        """
        return self.__class__(
                    [ -c for c in self.__coefficients ]
                )
    
    
    def __mul__(self, other):
        """
        Return the product of @p self and @p other. The infix operator @c *
        calls this method.
        
        Two polynomials are multiplied through convolution of their
        coefficients: for polynomials @f$ \sum_{k} a_{k}x^{k} @f$ and
        @f$ \sum_{k} b_{k}x^{k} @f$, their product is the polynomial
        @f$ \sum_{k} \sum_{j=0}^{k}(a_{j} + b_{k-j})x^{k} @f$.
        """
        # Initialize result as list of all zeros
        zero = self._coefficient_field.zero()
        # Add 2 because degrees count from 0.
        result = [ zero ] * (self.degree() + other.degree() + 2)
        
        for i, x in enumerate(self.__coefficients):
            for j, y in enumerate(other.__coefficients):
                result[i + j]  +=  x * y 
        
        return self.__class__( result )


    def __divmod__(self, other):
        """
        Return the quotient and remainder of @p self divided by @p other.
        
        The built-in method @c divmod() calls this method.  For the returned
        values @c quotient and @c remainder, we have
        @code
        self == quotient * other + remainder
        @endcode
        """
        # Lists will be modified, so copy them
        dividend = self.__coefficients[:]
        divisor = other.__coefficients[:]
        n = other.degree()
        
        zero = self._coefficient_field.zero()
        quotient = [ zero ] * (self.degree() - n + 1)
        
        for k in reversed(range( 0, len(quotient) )):
            quotient[k] = dividend[n + k] / divisor[n]
            for j in range(k, n + k):
                dividend[j] -= quotient[k] * divisor[j - k]
    
        remainder = dividend[ 0 : n ]
        
        return self.__class__( quotient ), \
                self.__class__( remainder )
    
    
    def __call__(self, point):
        """
        Return the polynomial function's value at @p point.
        
        This method evaluates the polynomial by replacing the indeterminate
        with @p point and summing the powers: for a polynomial
        @f$ p(x) = \sum_{k} a_{k}x^{k} @f$, the value at @f$ c @f$ is
        @f$ p(c) = \sum_{k} a_{k}c^{k} @f$.
        """
        return sum( [ c * point**i for i, c in enumerate(self.__coefficients) ] )
    
    
    #- Class Methods----------------------------------------------------------- 
    
    @classmethod
    def coefficient_field(cls):
        """
        Return the coefficient field.
        """
        return cls._coefficient_field


    @classmethod
    def zero(cls):
        """
        Return the polynomial ring's neutral element of addition: the zero
        polynomial.
        """
        return cls( cls._coefficient_field.zero() )
    
    
    @classmethod
    def one(cls):
        """
        Return the polynomial ring's neutral element of multiplication: the
        constant polynomial one.
        """
        return cls( cls._coefficient_field.one() )


    #- Auxiliary Functions ---------------------------------------------------- 

    def __remove_leading_zeros(self):
        """
        Remove all leading zeros from the list of coefficients.  This might
        empty the list completely. 
        """
        while len( self.__coefficients ) > 0 \
                and not self.__coefficients[-1]:
            self.__coefficients.pop()


    @staticmethod
    def __pad_and_zip(list1, list2, padding_element):
        """
        Combine @p list1 and @p list1 into a single list of pairs.  If the
        lists have different lengths, pad the shorter list with
        @p padding_element: fill in elements until they have the same length.
        """
        max_length = max( len(list1), len(list2) )
        padded_list1 = list1 + ( [padding_element] * (max_length - len(list1)) )
        padded_list2 = list2 + ( [padding_element] * (max_length - len(list2)) )
        return zip( padded_list1, padded_list2 ) class QuotientRing( CommutativeRing, metaclass=template( "_ring", "_modulus" ) ):
    """
    A ring of residue classes of ring elements modulo a fixed ring element; the
    residue classes (quotient ring elements) support infix notation for ring
    operations, as well as mixed argument types.
    
    This is a template class that must be instantiated with the source ring
    and modulus.  Use it, for example, as follows:
    @code
    # Instantiate the template; Z4 is a class (here: integer congruence classes)
    Z4 = QuotientRing( rings.integers.naive.Integers, 4 )
    x = Z4(1)      # Create a residue class: (1 mod 4)
    y = Z4(3)      # Another residue class: (2 mod 4)
    z = x+y        # z is (0 mod 4)
    z == Z4(2) * 2 # This true because 4 == 0 modulo 4
    type(z) is Z4  # This is also true
    @endcode
    
    Other examples of QuotientRing template instances are:
    @code
    Z17 = QuotientRing( rings.integers.naive.Integers, 17 )    # Z/17Z
    
    # The ring (Z/2Z[x]) / (x^3 + x + 1) is isomorphic to GF8,
    # the field with 8 elements
    GF2  = fields.finite.naive.FiniteField( 2 )
    GF2x = rings.polynomials.naive.Polynomials( GF2 )
    GF8 = QuotientRing( GF2x, GF2x(1, 1, 0, 1) ) 
    @endcode

    A quotient ring, factor ring, or residue class ring, is a ring
    @f$ R/(m) @f$, where @f$ R @f$ is the source ring, and @f$ (m) @f$
    denotes the ideal generated by the element @f$ m\in R @f$.  The ring
    operations carry over to quotient classes in a natural way: if
    @f$ x + y = z @f$ in @f$ R @f$, then @f$ [x] + [y] = [z] @f$
    in @f$ R/(m) @f$ (by the canonical homomorphism).

    With the quotient ring operations being defined in terms of source ring
    operations, elements of the source ring must implement:
    - __bool__(): Zero testing (@c True if not zero)
    - __eq__(): Equality testing with the @c == operator (@c True if equal)
    - __add__(): Addition with the @c + operator; @c self is the left summand
    - __neg__(): Negation with the @c - unary minus operator (the additive inverse)
    - __mul__(): Multiplication with the @c * operator; @c self is the
      left factor
    - __divmod__(): Division with remainder; @c self is the dividend (left element)
    
    @note  The implementation emphasizes simplicity over speed; it omits
           possible optimizations.
    
    @note  The class uses the operand_casting() decorator: @c other operands in
           binary operations will first be treated as QuotientRing elements.
           If that fails, the operation will be repeated after @p other was fed
           to the constructor __init__().  If that fails, too, then the
           operation returns @c NotImplemented (so that @p other.__rop__()
           might be invoked).
    
    @see   For example, Robinson, Derek J. S.,
           "An Introduction to Abstract Algebra", p. 106.
    """
    
    #- Instance Methods ----------------------------------------------------------- 
    
    def __init__(self, representative):
        """
        Construct a new residue class @p representative modulo modulus().
        
        If the @p representative already is an element of this QuotientRing
        class, then the new element is a copy of @p representative.
        """
        if isinstance( representative, self.__class__ ):
            self.__remainder = representative.__remainder
        elif isinstance( representative, self._modulus.__class__ ):
            self.__remainder = representative % self._modulus
        else:
            m = self._modulus
            self.__remainder = m.__class__( representative ) % m


    def remainder(self):
        """
        Return the remainder of the residue class (QuotientRing element)
        @p self.  This is an element of the source ring(), not a residue class.
        """
        return self.__remainder
    
    
    def __bool__(self):
        """
        Test whether the residue class (QuotientRing element) is non-zero:
        return @c True if, and only if, the remainder modulo modulus() is
        non-zero. Return @c False if the remainder is zero.
        
        Implicit conversions to boolean (truth) values use this method, for
        example when @c x is an element of a QuotientRing:
        @code
        if x:
            do_something()
        @endcode
        """
        return bool( self.__remainder )


    def __eq__(self, other):
        """
        Test whether another residue class (QuotientRing element) @p other is
        equivalent to @p self; return @c True if that is the case.  The infix
        operator @c == calls this method.
        
        Two residue classes @f$ [x], [y] @f$ are equal if, and only if, the
        difference of two representatives is a multiple of the modulus():
        @f$ x-y = m\cdot z @f$. 
        """
        # The representative is always the remainder; an equality test suffices
        return self.__remainder == other.remainder()


    def __add__(self, other):
        """
        Return the sum of @p self and @p other. The infix operator @c + calls
        this method.
        
        The sum of two residue classes (QuotientRing elements) @f$ [x], [y] @f$
        is the residue class @f$ [x + y] @f$. 
        """
        return self.__class__(
                        self.__remainder + other.remainder()
                    )
    

    def __neg__(self):
        """
        Return the additive inverse of @p self, which is @f$ [-x] @f$
        for a residue class (QuotientRing element) @f$ [x] @f$. The negation
        operator @c -x (unary minus) calls this method.
        """
        return self.__class__( -self.__remainder )


    def __mul__(self, other):
        """
        Return the product of @p self and @p other. The infix operator @c *
        calls this method.
        
        The product of two residue classes (QuotientRing elements)
        @f$ [x], [y] @f$ is the residue class @f$ [x \cdot y] @f$. 
        """
        return self.__class__(
                        self.__remainder * other.remainder()
                    )


    def __truediv__(self, other):
        """
        Return the quotient of @p self and @p other: multiply @p self with
        @c other.multiplicative_inverse(), so that
        @code
        (self / other) * other == self
        @endcode
        The infix operator @c / calls this method.

        @exception ZeroDivisionError   if @p other is not a unit, that is, has
                                       no multiplicative inverse.
        """
        return self * other.multiplicative_inverse()
        
    def __rtruediv__(self, other):
        """
        Return the quotient of @p other and @p self: multiply @p other with
        @c self.multiplicative_inverse(), so that
        @code
        (other / self) * self == other
        @endcode
        The infix operator @c / calls this method if other.__truediv__()
        returned @c NotImplemented.

        @exception ZeroDivisionError   if @p other is not a unit, that is, has
                                       no multiplicative inverse.
        """
        return other * self.multiplicative_inverse()
    

    def multiplicative_inverse(self):
        """
        Return an residue class (QuotientRing element) @c n such that
        @c n * self is one().
    
        @exception ZeroDivisionError   if @p self is not a unit, that is, has
                                       no multiplicative inverse.
        """
        if not self.__remainder:
            raise ZeroDivisionError

        inverse, ignore, gcd = \
            extended_euclidean_algorithm( self.__remainder, self._modulus )
    
        if gcd == self._ring.one():
            return self.__class__( inverse )
        else:
            message = "element has no inverse: representative and modulus " \
                      "are not relatively prime"
            raise ZeroDivisionError( message ) 


    #- Class Methods----------------------------------------------------------- 
    
    @classmethod
    def modulus(cls):
        """
        Return the quotient ring's modulus; this is an element of the
        source ring(), not a residue class (QuotientRing element).
        """
        return cls._modulus


    @classmethod
    def ring(cls):
        """
        Return the source ring.
        """
        return cls._ring
    
    
    @classmethod
    def zero(cls):
        """
        Return the quotient ring's neutral element of addition: the residue
        class (QuotientRing element) of ring().zero()
        """
        return cls( cls._ring.zero() )
    
    
    @classmethod
    def one(cls):
        """
        Return the quotient ring's neutral element of multiplication: the
        residue class (QuotientRing element) of ring().one()
        """
        return cls( cls._ring.one() ) def gcd(u, v):
    """
    Determine the greatest common divisor of @p u and @p v.
    
    @note      If @p u and @p v are Polynomials, then the returned polynomial
               is not necessarily monic; that means it may have a leading
               coefficient other than one. To check whether two polynomials are
               relatively prime, test whether the result has degree() zero.
    
    @see       extended_euclidean_algorithm() for a description of working
               input types.
    """
    return extended_euclidean_algorithm( u, v )[2] class Integers(int):
    """
    The ring of integers.  This wrapper class extends Python's built-in integer
    type @c int with the interface provided by the other ring classes.
    
    @see   rings.polynomials.naive.Polynomials,
           rings.quotients.naive.QuotientRing
    """
    
    @staticmethod
    def zero():
        """
        Return @c 0, the neutral element of integer addition.
        """
        return 0

    @staticmethod
    def one():
        """
        Return @c 1, the neutral element of integer multiplication.
        """
        return 1 class QuotientRing( CommutativeRing, metaclass=template( "_ring", "_modulus" ) ):
    """
    A ring of residue classes of ring elements modulo a fixed ring element; the
    residue classes (quotient ring elements) support infix notation for ring
    operations, as well as mixed argument types.
    
    This is a template class that must be instantiated with the source ring
    and modulus.  Use it, for example, as follows:
    @code
    # Instantiate the template; Z4 is a class (here: integer congruence classes)
    Z4 = QuotientRing( rings.integers.naive.Integers, 4 )
    x = Z4(1)      # Create a residue class: (1 mod 4)
    y = Z4(3)      # Another residue class: (2 mod 4)
    z = x+y        # z is (0 mod 4)
    z == Z4(2) * 2 # This true because 4 == 0 modulo 4
    type(z) is Z4  # This is also true
    @endcode
    
    Other examples of QuotientRing template instances are:
    @code
    Z17 = QuotientRing( rings.integers.naive.Integers, 17 )    # Z/17Z
    
    # The ring (Z/2Z[x]) / (x^3 + x + 1) is isomorphic to GF8,
    # the field with 8 elements
    GF2  = fields.finite.naive.FiniteField( 2 )
    GF2x = rings.polynomials.naive.Polynomials( GF2 )
    GF8 = QuotientRing( GF2x, GF2x(1, 1, 0, 1) ) 
    @endcode

    A quotient ring, factor ring, or residue class ring, is a ring
    @f$ R/(m) @f$, where @f$ R @f$ is the source ring, and @f$ (m) @f$
    denotes the ideal generated by the element @f$ m\in R @f$.  The ring
    operations carry over to quotient classes in a natural way: if
    @f$ x + y = z @f$ in @f$ R @f$, then @f$ [x] + [y] = [z] @f$
    in @f$ R/(m) @f$ (by the canonical homomorphism).

    With the quotient ring operations being defined in terms of source ring
    operations, elements of the source ring must implement:
    - __bool__(): Zero testing (@c True if not zero)
    - __eq__(): Equality testing with the @c == operator (@c True if equal)
    - __add__(): Addition with the @c + operator; @c self is the left summand
    - __neg__(): Negation with the @c - unary minus operator (the additive inverse)
    - __mul__(): Multiplication with the @c * operator; @c self is the
      left factor
    - __divmod__(): Division with remainder; @c self is the dividend (left element)
    
    @note  The implementation emphasizes simplicity over speed; it omits
           possible optimizations.
    
    @note  The class uses the operand_casting() decorator: @c other operands in
           binary operations will first be treated as QuotientRing elements.
           If that fails, the operation will be repeated after @p other was fed
           to the constructor __init__().  If that fails, too, then the
           operation returns @c NotImplemented (so that @p other.__rop__()
           might be invoked).
    
    @see   For example, Robinson, Derek J. S.,
           "An Introduction to Abstract Algebra", p. 106.
    """
    
    #- Instance Methods ----------------------------------------------------------- 
    
    def __init__(self, representative):
        """
        Construct a new residue class @p representative modulo modulus().
        
        If the @p representative already is an element of this QuotientRing
        class, then the new element is a copy of @p representative.
        """
        if isinstance( representative, self.__class__ ):
            self.__remainder = representative.__remainder
        elif isinstance( representative, self._modulus.__class__ ):
            self.__remainder = representative % self._modulus
        else:
            m = self._modulus
            self.__remainder = m.__class__( representative ) % m


    def remainder(self):
        """
        Return the remainder of the residue class (QuotientRing element)
        @p self.  This is an element of the source ring(), not a residue class.
        """
        return self.__remainder
    
    
    def __bool__(self):
        """
        Test whether the residue class (QuotientRing element) is non-zero:
        return @c True if, and only if, the remainder modulo modulus() is
        non-zero. Return @c False if the remainder is zero.
        
        Implicit conversions to boolean (truth) values use this method, for
        example when @c x is an element of a QuotientRing:
        @code
        if x:
            do_something()
        @endcode
        """
        return bool( self.__remainder )


    def __eq__(self, other):
        """
        Test whether another residue class (QuotientRing element) @p other is
        equivalent to @p self; return @c True if that is the case.  The infix
        operator @c == calls this method.
        
        Two residue classes @f$ [x], [y] @f$ are equal if, and only if, the
        difference of two representatives is a multiple of the modulus():
        @f$ x-y = m\cdot z @f$. 
        """
        # The representative is always the remainder; an equality test suffices
        return self.__remainder == other.remainder()


    def __add__(self, other):
        """
        Return the sum of @p self and @p other. The infix operator @c + calls
        this method.
        
        The sum of two residue classes (QuotientRing elements) @f$ [x], [y] @f$
        is the residue class @f$ [x + y] @f$. 
        """
        return self.__class__(
                        self.__remainder + other.remainder()
                    )
    

    def __neg__(self):
        """
        Return the additive inverse of @p self, which is @f$ [-x] @f$
        for a residue class (QuotientRing element) @f$ [x] @f$. The negation
        operator @c -x (unary minus) calls this method.
        """
        return self.__class__( -self.__remainder )


    def __mul__(self, other):
        """
        Return the product of @p self and @p other. The infix operator @c *
        calls this method.
        
        The product of two residue classes (QuotientRing elements)
        @f$ [x], [y] @f$ is the residue class @f$ [x \cdot y] @f$. 
        """
        return self.__class__(
                        self.__remainder * other.remainder()
                    )


    def __truediv__(self, other):
        """
        Return the quotient of @p self and @p other: multiply @p self with
        @c other.multiplicative_inverse(), so that
        @code
        (self / other) * other == self
        @endcode
        The infix operator @c / calls this method.

        @exception ZeroDivisionError   if @p other is not a unit, that is, has
                                       no multiplicative inverse.
        """
        return self * other.multiplicative_inverse()
        
    def __rtruediv__(self, other):
        """
        Return the quotient of @p other and @p self: multiply @p other with
        @c self.multiplicative_inverse(), so that
        @code
        (other / self) * self == other
        @endcode
        The infix operator @c / calls this method if other.__truediv__()
        returned @c NotImplemented.

        @exception ZeroDivisionError   if @p other is not a unit, that is, has
                                       no multiplicative inverse.
        """
        return other * self.multiplicative_inverse()
    

    def multiplicative_inverse(self):
        """
        Return an residue class (QuotientRing element) @c n such that
        @c n * self is one().
    
        @exception ZeroDivisionError   if @p self is not a unit, that is, has
                                       no multiplicative inverse.
        """
        if not self.__remainder:
            raise ZeroDivisionError

        inverse, ignore, gcd = \
            extended_euclidean_algorithm( self.__remainder, self._modulus )
    
        if gcd == self._ring.one():
            return self.__class__( inverse )
        else:
            message = "element has no inverse: representative and modulus " \
                      "are not relatively prime"
            raise ZeroDivisionError( message ) 


    #- Class Methods----------------------------------------------------------- 
    
    @classmethod
    def modulus(cls):
        """
        Return the quotient ring's modulus; this is an element of the
        source ring(), not a residue class (QuotientRing element).
        """
        return cls._modulus


    @classmethod
    def ring(cls):
        """
        Return the source ring.
        """
        return cls._ring
    
    
    @classmethod
    def zero(cls):
        """
        Return the quotient ring's neutral element of addition: the residue
        class (QuotientRing element) of ring().zero()
        """
        return cls( cls._ring.zero() )
    
    
    @classmethod
    def one(cls):
        """
        Return the quotient ring's neutral element of multiplication: the
        residue class (QuotientRing element) of ring().one()
        """
        return cls( cls._ring.one() ) class FiniteField(QuotientRing, _ring=Integers):
    """
    The finite field with @f$ p^k @f$ elements; the field operations on the
    element instances support the natural infix syntax and allow mixing in
    integer arguments.
    
    This is a template class that must be instantiated with the field size.
    For example:
    @code
    # Instantiate the template; GF23 is a class.
    GF23 = FiniteField(23)
    x = GF23(3)    # Create a field element: x is the quotient class (3 mod 23)
    y = GF23(-2)   # Equivalently: y = GF23(21) because 21 == -2  (mod 23)
    z = x**2 + 2*y - x/y   # z == 15 (mod 23); note that -2*11 == 1 (mod 23)
    type(x) is GF23        # This is True
    @endcode
    
    The implementation is held as simple as possible. Its focus lies on ease
    of understanding, not on performance. Thus it is badly suited for larger
    computations. The template inherits from
    @c rings.quotients.naive.QuotientRing and specializes the @c _ring template
    parameter to @c rings.integers.naive.Integers 
    
    @note      The field size must be a prime. The template will, however,
               accept compound numbers without hesitation. In this case, the
               class behaves like the quotient ring
               @f$ \mathbb{Z}/n\mathbb{Z} @f$.
    
    @see       rings.quotients.naive.QuotientRing
    
    @author    Peter Dinges <pdinges@acm.org>
    """

    @classmethod
    def characteristic(cls):
        """
        Return the field characteristic @f$ p @f$.
        """
        return cls._modulus
    
    @classmethod
    def power(cls):
        """
        Return the power @f$ k @f$ of the characteristic that yields the field
        size: @f$ q = p^k @f$, where @f$ p @f$ is the field characteristic.
        
        @note  In the current implementation, @f$ k @f$ is always 1.
        """
        return 1
    
    @classmethod
    def size(cls):
        """
        Return the number of elements in the field.
        """
        return cls.characteristic() ** cls.power()

    @classmethod
    def elements(cls):
        """
        Return a list of all field elements in ascending order
        @f$ (0, 1, 2, ..., p-1) @f$ (mod p).
        
        @note  The method populates the complete list, so this operation might
               be expensive for large field characteristics.
        """
        return [ cls(i) for i in range(0, cls.size()) ] class EllipticCurve( metaclass=template("_field", "_A", "_B") ):
    """
    An elliptic curve; the points support additive infix notation for group
    operations, as well as multiplication by integers.
    
    This is a template class that must be instantiated with the field and the
    parameters @c A and @c B.  Use it, for example, as follows:
    @code
    # Instantiate the elliptic curve template as y^2 = x^3 + 3x + 4 over GF(7)
    E = EllipticCurve( FiniteField(7), 3, 4 )

    # Create two points and add them
    P = E(1, 1)
    Q = E(2, 5)
    
    P + Q == E(6, 0)   # Point addition
    3*P == E(0, 5)     # Multiplication by integers
    
    O = PointAtInfinity()
    P == P + O         # The point at infinity is the neutral element
    type(P) is E       # This is also true
    @endcode
    
    An elliptic curve is the set of all points @f$ E = \bigl\{ (x, y) \mid y^2
    = x^3 + Ax + B \bigr\} \cup \{\mathcal{O}\} @f$, where @f$ \mathcal{O} @f$
    is a special element, the point at infinity.  The points form an abelian
    group with the point at infinity as neutral element.
    
    The relation @f$ y^2 = x^3 + Ax + B @f$ between the (Cartesian) coordinates
    of the points @f$ (x, y) \in E @f$ is called the fundamental relation.

    @see   Charlap, Leonard S. and Robbins, David P., "CRD Expositroy Report 31:
           An Elementary Introduction to Elliptic Curves", 1988, chapter 1
    """
    
    #- Instance Methods ----------------------------------------------------------- 
    
    def __init__(self, x, y):
        """
        Construct a new point on the elliptic curve with coordinates @p x
        and @p y.
        
        The coordinates must be connected by the fundamental relation
        @f$ y^2 = x^3 + Ax + B @f$.
        
        @param x   The first coordinate of the point.  It will be interpreted
                   as an element of field(), which might result in a ValueError
                   or TypeError.
        @param y   The second coordinate of the point.  It will be interpreted
                   as an element of field(), which might result in a ValueError
                   or TypeError.
        
        @exception AssertionError  if the fundamental relation does not hold.
        @exception ValueError      if an argument cannot be cast as element
                                   of field(). 
        @exception TypeError       same as @c ValueError. 
        """
        self.__x = self._field( x )
        self.__y = self._field( y )

        A, B = self.parameters()
        x, y = self.__x, self.__y
        assert y ** 2 == x ** 3  +  A * x  +  B, \
            "point ({x}, {y}) is not on the curve".format(x=x, y=y)
    
    
    def x(self):
        """
        Return the first (Cartesian) coordinate of the point @f$ (x,y) @f$
        on the curve.
        
        @note  The fundamental relation links the coordinates:
               @f$ y^2 = x^3 + Ax + B @f$.
        """
        return self.__x
    
    
    def y(self):
        """
        Return the second (Cartesian) coordinate of the point @f$ (x,y) @f$
        on the curve.
        
        @note  The fundamental relation links the coordinates:
               @f$ y^2 = x^3 + Ax + B @f$.
        """
        return self.__y
    
    
    def is_infinite(self):
        """
        Test whether the point is infinite or not: always return @c False for
        the point is finite.  The point at infinity returns @c True.
        """
        return False

    
    def __bool__(self):
        """
        Test whether the point is infinite or not: always return @c True, for
        the point is finite.  The point at infinity returns @c False.
        
        Implicit conversions to boolean (truth) values use this method, for
        example when @c P is an element of EllipticCurve:
        @code
        if P:
            do_something()
        @endcode
        """
        return True
    
        
    def __eq__(self, other):
        """
        Test whether another point @p other is equal to @p self; return
        @c True if that is the case.  The infix operator @c == calls
        this method.
        
        Two points are equal if, and only if, both of their (Cartesian)
        coordinates are equal or if they both are infinite.
        """
        if other.is_infinite():
            return False
        else:
            return self.__x == other.x() and self.__y == other.y()
    
    
    def __neq__(self, other):
        """
        Test whether another point @p other is different from @p self; return
        @c True if that is the case.  The infix operator @c != calls
        this method.
        
        Two points are different if, and only if, they differ in at least one
        of their (Cartesian) coordinates. 
        """
        return not self == other
    
    
    def __add__(self, other):
        """
        Return the sum of @p self and @p other.  The infix operator @c + calls
        this method.
        
        The geometric interpretation of point addition is to draw a line between
        the points and mirror the point at the third intersection of line and
        curve on the y-axis.  However, the justification why this works is
        rather complicated.  Please see some introductory text for the formulas
        and the respective proofs.
        
        @note  Addition on elliptic curves is commutative.  The neutral element
               is the point at infinity (compare PointAtInfinity).
        
        @see   For example Washington, Lawrence C. "Elliptic Curves:
               Number Theory and Cryptography", second edition, CRC Press 2008,
               chapter 2.
        """
        if other.is_infinite():
            return self
        elif self == -other:
            return PointAtInfinity()
        else:
            if self.__x == other.x():
                return self.__double__()
            else:
                return self.__generic_add__( other )

                
    def __generic_add__(self, other):
        """
        Generic addition of points @p self and @p other: the points are neither
        identical nor inverses.
        
        @note  This method should not be called directly.  Instead use the
               infix operator @c +, which calls __add__(). 
        
        @note  Putting the cases in separate methods makes call profiles more
               expressive.
        
        @see __add__()
        """
        gamma = (other.y() - self.__y) / (other.x() - self.__x)
        u = -self.__x - other.x() +  gamma ** 2
        v = -self.__y - gamma * (u - self.__x)
        
        return self.__class__(u, v)


    def __double__(self):
        """
        Point doubling: add a point to itself.
        
        @note  This method should not be called directly.  Instead use the
               infix operator @c +, which calls __add__().
        
        @note  Putting the cases in separate methods makes call profiles more
               expressive.
        
        @see __add__()
        """
        A, B = self.parameters()
        delta = (3 * self.__x ** 2  + A) / (2 * self.__y)
        u = -self.__x - self.__x +  delta ** 2
        v = -self.__y - delta * (u - self.__x)
        
        return self.__class__(u, v)
    
    
    def __neg__(self):
        """
        Return the additive inverse (the negative) of @p self.
        
        The additive inverse of a point @f$ (x,y) @f$ on the curve is
        @f$ (x, -y) @f$.
        """
        return self.__class__(self.__x, -self.__y)
    
    
    def __sub__(self, other):
        """
        Return the difference between @p self and @p other.  The infix operator
        @c - calls this method.
        
        This is the same as @p self + (-other).
        
        @see   __add__()
        """
        return self + (-other)
    
    
    def __mul__(self, n):
        """
        Multiplication with integers: adding @p n copies of the point @p self;
        the infix operator @c * calls this method.
        
        This method is used when self was the left factor.

        @param n   A number object that can be interpreted as an integer.  It
                   determines how often the point will be added to itself.
        
        @exception ValueError  if @p n cannot be cast to @c int().
        @exception TypeError   same as @c ValueError.  
        """
        n = int(n)
        if n == 0:
            return PointAtInfinity()
        
        point = self
        for i in range(1, n):
            point += self
        
        return point
    
    def __rmul__(self, other):
        """
        Multiplication with integers: adding @p n copies of the point @p self;
        the infix operator @c * calls this method.
        
        This method is used when self was the right factor.

        @param n   A number object that can be interpreted as an integer.  It
                   determines how often the point will be added to itself.
        
        @exception ValueError  if @p n cannot be cast to @c int().
        @exception TypeError   same as @c ValueError.  
        """
        # Multiplication with integers is always commutative.
        return self * other
    
    
    #- Class Methods----------------------------------------------------------- 
    
    @classmethod
    def field(cls):
        """
        Return the field over which the curve was defined.
        """
        return cls._field
    
    
    @classmethod
    def parameters(cls):
        """
        Return the parameter pair @f$ (A, B) @f$ of the curve.
        
        The parameters determine the fundamental relation
        @f$ y^2 = x^3 + Ax + B @f$.
        """
        return (cls._A, cls._B)
    
    
    @classmethod
    def is_singular(cls):
        """
        Test whether the curve is singular or not: return @c True if, and
        only if, the curve is singular.
        
        A curve is singular if the polynomial @f$ x^3 + Ax + B @f$ has
        repeated roots.
        """
        # A curve is singular if and only if its determinant is 0
        # (in fields of characteristic neither 2 nor 3).
        return not 4 * cls._A**3  + 27 * cls._B ** 2 class AlgorithmRunner:
    """
    An AlgorithmRunner allows fast creation of a feature-rich command line
    program that applies an algorithm to sets of parameters while recording
    the execution profile.
    
    For example, say you have an algorithm (that is, a callable) @c f(x,y) and
    want to apply it to every parameter set listed in the file @c params.lst.
    Furthermore suppose you want to record the calling profile (see the
    @c cProfile module), the execution time, and the result of each
    application. Then a script containing the following code is your program
    (@c prog.py):
    @code
    def f(x, y, output):
        # Do calculations here
        print( "The ultimate truth is:", result, file=output )
        return result
    
    r = AlgorithmRunner( f )
    r.run()
    @endcode
    You then invoke
    @code
    python prog.py -d profiles/ -p -i params.lst
    @endcode
    and it will apply @c f() to each line of @c params.lst and store the
    respective profile and timing information in the @c profiles/ directory.
    Use the "-h" command line switch to get a list of all available options.
    
    AlgorithmRunner uses ParallelParser to read the input files. It therefore
    supports multiple program instances working on the same parameter list.
    Each line will be read by only one instance. 
    
    Send the signal @c SIGUSR1 to an AlgorithmRunner instance to "soft" kill
    it. It will then properly shut down, unregistering its ParallelParser etc.,
    all of which does not happen if it is just killed.
    
    @see ParallelParser and DataRecorder
    """
    
    def __init__(self, algorithm, arguments=sys.argv[1:], algorithm_version="<unknown>" ):
        """
        Construct a new AlgorithmRunner for executing the given @p algorithm.
        
        @param     algorithm   A callable that takes @f$ n+1 @f$ arguments,
                               where @f$ n @f$ is the number of entries per line
                               in parameter lists (or given on the command line).
                               The last argument is @c output, which gives the
                               algorithm a way to write logging information.
                               The return value of @p algorithm will be stored
                               alongside the timing information.
        @param     arguments   A list of command line argument strings for the
                               runner.
        @param     algorithm_version   The algorithm version will be stored
                                       together with the timing information.
                                       Use it to recall which revision of an
                                       algorithm was used in a particular
                                       profile.
        """
        self.__algorithm = algorithm
        self.__algorithm_version = algorithm_version
        options, arguments = self._parse_arguments( arguments, algorithm_version )
        
        # __input is a list of pairs (<name>, <iterable>);
        # <iterable> is expected to return pairs (<item_number>, <item>).
        # See run().
        self.__input = []
        if arguments:
            self.__input.append( ( "<stdin>", [ (0, tuple( arguments ) ) ] ) )
        # Fail early: immediately try to open the file
        if options.input_file:
            input_parser = ParallelParser( options.input_file )
            self.__input.append( ( options.input_file, input_parser ) )
        
        # Initialize the remaining attributes.
        self._open_output( options.output_file )

        self.__create_profile = options.create_profile
        self.__profile_directory = options.profile_directory
        self.__timelimit = options.timelimit
        
        if self.__create_profile:
            self._ensure_directory( options.profile_directory )


    def run(self):
        """
        Start execution and apply the algorithm to all provided parameters.
        """
        if not self.__input:
            self.__parser.print_usage()
        
        # Listen for signals to allow the user to interrupt the execution.
        signal.signal( signal.SIGUSR1, self._raise_signal )
        
        terminated = False    
        for name, iterable in self.__input:
            for item_number, item in iterable:
                try:
                    # Raise a TimeOutException if the pass takes too long
                    if self.__timelimit > 0:
                        signal.signal( signal.SIGALRM, self._raise_signal )
                        signal.alarm( self.__timelimit )
                        
                    if self.__create_profile:
                        dump_filename = self._generate_dump_name( *item )
                        recorder = DataRecorder(
                                        dump_filename,
                                        self.__profile_directory,
                                        
                                    )
                    try:
                        result = self.__algorithm( *item, output=self.__output)
                    
                    except TimeOutException:
                        terminated = "timeout after {0}s".format( self.__timelimit )
                    
                    except SoftKillException:
                        terminated = "killed by user signal"
                    
                    if self.__create_profile:
                        recorder.stop()
                        extra_information = {
                                     "algorithm": self.__algorithm.__name__,
                                     "version": self.__algorithm_version,
                                     "input": item,
                                 }

                        if terminated:
                            extra_information["terminated"] = terminated
                        else:
                            extra_information["result"] = result

                        recorder.dump_data( extra_information )
    
                except TypeError:
                    message = ">>> Algorithm input in '{name}':{i} has the " \
                              "wrong type. Ignoring it."
                    print( message.format(
                                  name = name,
                                  i = item_number + 1
                              ),
                              file = sys.stderr
                          )
                except IOError:
                    pass # output failure
                
                if terminated:
                    message = ">>> Terminated prematurely. Reason: {0}"
                    print( message.format( terminated ), file=sys.stderr ) 
                    return 1
        
        # Stop listening for signals.
        signal.signal( signal.SIGUSR1, signal.SIG_DFL )
        signal.signal( signal.SIGALRM, 0 )


    def _ensure_directory(self, path):
        """
        Make sure that the directory @p path exists.
        """
        if not os.path.exists( path ):
            os.makedirs( path )
            
        elif not os.path.isdir( path ):
            message = "file already exists and is not a path: '{0}'"
            raise IOError( message.format( path ) )


    # TODO: Distinguish results, errors, and log
    def _open_output(self, filename, default = sys.stdout):
        """
        Make sure that @p filename is ready for writing.
        """
        if filename and os.path.exists( filename ):
            raise IOError( "output file already exists" )
        
        if filename:
            self.__output = open( filename, mode="wt" )
        else:
            self.__output = default

    
    def _generate_dump_name(self, *input_item):
        """
        Return a generic file name for profiling information.
        """
        name = "{date}-{args}".format(
                   date = datetime.now().strftime("%Y%m%d_%H%M%S"),
                   args = "_".join( input_item )
               )
        return name


    def _parse_arguments( self, arguments, algorithm_version ):
        """
        Parse the command line arguments and set the options accordingly.
        """
        usage_string = """%prog <algorithm input>"""
        version_string ="""%prog {0}""".format( algorithm_version )
        
        parser = OptionParser( usage = usage_string, version = version_string )
        
        parser.add_option( "-t",
                           "--timelimit",
                           dest="timelimit",
                           type="int",
                           default=0,
                           metavar="SECONDS",
                           help="stop execution if an algorithm pass takes "
                           "longer than SECONDS seconds."
                       )

        io_group = OptionGroup(parser, "Input and Output")
        io_group.add_option( "-i", "--input-file", metavar="FILE",
                             dest="input_file",
                             type="string",
                             default=None,
                             help="read input line by line from FILE")
        
        io_group.add_option( "-o", "--output-file", metavar="FILE",
                             dest="output_file",
                             type="string",
                             default=None,
                             help="write results to FILE")

        # TODO: Add overwrite option.
        parser.add_option_group( io_group )

        
        profiling_group = OptionGroup(parser, "Profiling")
        profiling_group.add_option( "-p", "--create-profile",
                           action="store_true", dest="create_profile",
                           default=False,
                           help="create a profile for each algorithm run"
                           )
        
        profiling_group.add_option( "-d", "--profile-directory", metavar="DIR",
                           dest="profile_directory",
                           type="string",
                           default=os.getcwd(),
                           help="store profiles in directory DIR")
        
        parser.add_option_group( profiling_group )
        
        self.__parser = parser
        return parser.parse_args( arguments )
    
    
    def _raise_signal( self, signum, frame ):
        """
        Asynchronous signal handler that raises the respective exceptions.
        
        @see TimeOutException, SoftKillException
        """
        if signum == signal.SIGALRM:
            raise TimeOutException
        elif signum == signal.SIGUSR1:
            raise SoftKillException
from elliptic_curves.l_torsion_group.naive import LTorsionGroup
from support.primes import inverse_primorial, primes_range
from support.quotients import solve_congruence_equations, representative_in_range
from rings.integers.naive import Integers
from rings.polynomials.naive import Polynomials
from rings.quotients.naive import QuotientRing
from support.rings import gcd
from rings.integers.naive import Integers
from rings.quotients.naive import QuotientRing
from fields.finite.naive import FiniteField
from elliptic_curves.naive import EllipticCurve
from support.running import AlgorithmRunner
import sys
# -*- coding: utf-8 -*-
# Copyright (c) 2010--2012  Peter Dinges <pdinges@acm.org>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.


"""
A naive implementation of Schoof's algorithm for counting the points of
elliptic curves over finite fields.

The implementation is meant for instructional purposes, not for actual counting
of, say, curves for cryptographic settings. It uses naive arithmetic and
emphasizes ease of understanding over performance; all speed-ups (such as
symmetry of search ranges, Hasse's theorem, etc.) are excluded.

@see Schoof, R.
"Elliptic Curves Over Finite Fields and the Computation of Square Roots mod p"
in Mathematics of Computation, Vol. 44, No. 170 (Apr. 1985), pp. 483--494.

@package   schoof.naive
@author    Peter Dinges <pdinges@acm.org>
"""


def frobenius_trace(curve):
    """
    Compute the trace of the Frobenius endomorphism for the given EllpiticCurve
    @p curve.
    
    This is an implementation of Schoof's original algorithm for counting the
    points of an elliptic curve over a finite field.
    
    @return    The trace @f$ t @f$ of the Frobenius endomorphism. The number of
               points on the curve then is @f$ q + 1 - t @f$, where @f$ q @f$
               is the size of the finite field over which the curve was defined.
    """
    # Initialize variables and parameters
    trace_congruences = []
    search_range = possible_frobenius_trace_range( curve.field() )
    upper_prime_bound = inverse_primorial(
                            len(search_range),
                            shunned = curve.field().characteristic()
                          )
    
    # Collect the congruence equations (avoid multivariate
    # polynomial arithmetic by handling 2-torsion separately)
    trace_congruences.append( frobenius_trace_mod_2( curve ) )

    torsion_group = LTorsionGroup( curve )

    for prime in primes_range( 3, upper_prime_bound+1 ):