def __init__(self, p, check=True, modulus=None):
        """
        Return a new finite field of order `p` where `p` is prime.

        INPUT:

        - ``p`` -- an integer at least 2

        - ``check`` -- bool (default: ``True``); if ``False``, do not
          check ``p`` for primality

        EXAMPLES::

            sage: F = FiniteField(3); F
            Finite Field of size 3
        """
        p = Integer(p)
        if check and not p.is_prime():
            raise ArithmeticError("p must be prime")
        self.__char = p
        # FiniteField_generic does nothing more than IntegerModRing_generic, and
        # it saves a non trivial overhead
        integer_mod_ring.IntegerModRing_generic.__init__(self, p, category=_FiniteFields)

        # If modulus is None, it will be created on demand as x-1
        # by the modulus() method.
        if modulus is not None:
            self._modulus = modulus

    def __getitem__(self,level):
        """
        Return the modular polynomial of given level, or an error if
        there is no such polynomial in the database.
        
        EXAMPLES:
            sage: DBMP = ClassicalModularPolynomialDatabase()  #optional requires database_hohel
            sage: f = DBMP[29]                                 #optional
            sage: f.degree()                                   #optional
            58
            sage: f.coefficient([28,28])                       #optional
            400152899204646997840260839128

            sage: DBMP[50]                                     #optional
            Traceback (most recent call last):
            ...
            RuntimeError: No database entry for modular polynomial of level 50
        """
        if self.model in ("Atk","Eta"):
            level = Integer(level)
            if not level.is_prime():
                raise TypeError, "Argument level (= %s) must be prime."%level 
        elif self.model in ("AtkCrr","EtaCrr"):
            N = Integer(level[0])
            if not N in (2,3,5,7,13):
                raise TypeError, "Argument level (= %s) must be prime."%N
        modpol = self._dbpath(level)
        try:
            coeff_list = _dbz_to_integer_list(modpol)
        except RuntimeError, msg:
            print msg
            raise RuntimeError, \
                  "No database entry for modular polynomial of level %s"%level

def hurwitz_kronecker_class_no_x( n, D):
    r"""
    OUTPUT
        The generalized Hurwitz-Kronecker class number $H_n(D)$
        as defined in [S-Z].

    INPUT
        n -- an integer $\ge 1$
        D -- an integer $\le 0$
    """
    n = Integer(n)
    D = Integer(D)
    if n <= 0:
        raise ValueError, "%s: must be an integer >=1" % n
    if D > 0:
        raise ValueError, "%s: must be an integer <= 0" % D
    g = D.gcd( n)
    a,b = _ab(g)
    h = g*b
    if 0 != D%h:
        return 0
    Dp = D//h
    H1 = Rational( gp.qfbhclassno( -Dp))
    if 0 == H1:
        return H1
    return g * Dp.kronecker( n//g) * H1

    def _validate(self, n):
        """
        Verify that n is positive and has at most 4095 digits.

        INPUT:
            n

        This function raises a ValueError if the two conditions listed above
        are not both satisfied.  It is here because GMP-ECM silently ignores
        all digits of input after the 4095th!

        EXAMPLES:
            sage: ecm = ECM()
            sage: ecm._validate(0)
            Traceback (most recent call last):
            ...
            ValueError: n must be positive
            sage: ecm._validate(10^5000)
            Traceback (most recent call last):
            ...
            ValueError: n must have at most 4095 digits
        """
        n = Integer(n)
        if n <= 0:
            raise ValueError, "n must be positive"
        if n.ndigits() > 4095:
            raise ValueError, "n must have at most 4095 digits"

    def nth_root(self, n):
        """
        Return an `n`-th root of ``self``.

        EXAMPLES::

            sage: F = GF(5).algebraic_closure()
            sage: t = F.gen(2) + 1
            sage: s = t.nth_root(15); s
            4*z6^5 + 3*z6^4 + 2*z6^3 + 2*z6^2 + 4
            sage: s**15 == t
            True

        .. TODO::

            This function could probably be made faster.

        """
        from sage.rings.integer import Integer
        F = self.parent()
        x = self._value
        n = Integer(n)
        l = self._level
        # In order to be smart we look for the smallest subfield that
        # actually contains the root.
        for d in n.divisors():
            xx = F.inclusion(l, d*l)(x)
            try:
                y = xx.nth_root(n, extend=False)
            except ValueError:
                continue
            return self.__class__(F, y)

        raise AssertionError('cannot find n-th root in algebraic closure of finite field')

    def iter_indefinite_forms(self) :
        fm = Integer(4 * self.__m)
        
        if self.__reduced :
            if self.__weak_forms :
                msq = self.__m**2
                for n in xrange(0, min(self.__m // 4 + 1, self.__bound)) :
                    for r in xrange( isqrt(fm * n - 1) + 1 if n != 0 else 0,
                                     isqrt(fm * n + msq + 1) ) :
                        yield (n, r)
            else :

                for r in xrange(0, min(self.__m + 1,
                                       isqrt((self.__bound - 1) * fm) + 1) ) :
                    if fm.divides(r**2) :
                        yield (r**2 // fm, r)
        else :
            if self.__weak_forms :
                msq = self.__m**2
                for n in xrange(0, self.__bound) :
                    for r in xrange( isqrt(fm * n - 1) + 1 if n != 0 else 0,
                                     isqrt(fm * n + msq + 1) ) :
                        yield (n, r)
            else :
                for n in xrange(0, self.__bound) :
                    if (fm * n).is_square() :
                        yield(n, isqrt(fm * n))

        raise StopIteration

    def __init__(self, q, name="a", modulus=None, repr="poly", cache=False):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: k.<a> = GF(2^3)
            sage: j.<b> = GF(3^4)
            sage: k == j
            False

            sage: GF(2^3,'a') == copy(GF(2^3,'a'))
            True
            sage: TestSuite(GF(2^3, 'a')).run()
        """
        self._kwargs = {}

        if repr not in ["int", "log", "poly"]:
            raise ValueError, "Unknown representation %s" % repr

        q = Integer(q)
        if q < 2:
            raise ValueError, "q  must be a prime power"
        F = q.factor()
        if len(F) > 1:
            raise ValueError, "q must be a prime power"
        p = F[0][0]
        k = F[0][1]

        if q >= 1 << 16:
            raise ValueError, "q must be < 2^16"

        import constructor

        FiniteField.__init__(self, constructor.FiniteField(p), name, normalize=False)

        self._kwargs["repr"] = repr
        self._kwargs["cache"] = cache

        self._is_conway = False
        if modulus is None or modulus == "conway":
            if k == 1:
                modulus = "random"  # this will use the gfq_factory_pk function.
            elif ConwayPolynomials().has_polynomial(p, k):
                from sage.rings.finite_rings.constructor import conway_polynomial

                modulus = conway_polynomial(p, k)
                self._is_conway = True
            elif modulus is None:
                modulus = "random"
            else:
                raise ValueError, "Conway polynomial not found"

        from sage.rings.polynomial.all import is_Polynomial

        if is_Polynomial(modulus):
            modulus = modulus.list()

        self._cache = Cache_givaro(self, p, k, modulus, repr, cache)

def QuadraticResidueCodeOddPair(n, F):
    """
    Quadratic residue codes of a given odd prime length and base ring
    either don't exist at all or occur as 4-tuples - a pair of
    "odd-like" codes and a pair of "even-like" codes. If n 2 is prime
    then (Theorem 6.6.2 in [HP]_) a QR code exists over GF(q) iff q is a
    quadratic residue mod n.

    They are constructed as "odd-like" duadic codes associated the
    splitting (Q,N) mod n, where Q is the set of non-zero quadratic
    residues and N is the non-residues.

    EXAMPLES::

        sage: codes.QuadraticResidueCodeOddPair(17,GF(13))
        (Linear code of length 17, dimension 9 over Finite Field of size 13,
         Linear code of length 17, dimension 9 over Finite Field of size 13)
        sage: codes.QuadraticResidueCodeOddPair(17,GF(2))
        (Linear code of length 17, dimension 9 over Finite Field of size 2,
         Linear code of length 17, dimension 9 over Finite Field of size 2)
        sage: codes.QuadraticResidueCodeOddPair(13,GF(9,"z"))
        (Linear code of length 13, dimension 7 over Finite Field in z of size 3^2,
         Linear code of length 13, dimension 7 over Finite Field in z of size 3^2)
        sage: C1 = codes.QuadraticResidueCodeOddPair(17,GF(2))[1]
        sage: C1x = C1.extended_code()
        sage: C2 = codes.QuadraticResidueCodeOddPair(17,GF(2))[0]
        sage: C2x = C2.extended_code()
        sage: C2x.spectrum(); C1x.spectrum()
        [1, 0, 0, 0, 0, 0, 102, 0, 153, 0, 153, 0, 102, 0, 0, 0, 0, 0, 1]
        [1, 0, 0, 0, 0, 0, 102, 0, 153, 0, 153, 0, 102, 0, 0, 0, 0, 0, 1]
        sage: C3 = codes.QuadraticResidueCodeOddPair(7,GF(2))[0]
        sage: C3x = C3.extended_code()
        sage: C3x.spectrum()
        [1, 0, 0, 0, 14, 0, 0, 0, 1]

    This is consistent with Theorem 6.6.14 in [HP]_.

    TESTS::

        sage: codes.QuadraticResidueCodeOddPair(9,GF(2))
        Traceback (most recent call last):
        ...
        ValueError: the argument n must be an odd prime
    """
    from sage.arith.srange import srange
    from sage.categories.finite_fields import FiniteFields

    if F not in FiniteFields():
        raise ValueError("the argument F must be a finite field")
    q = F.order()
    n = Integer(n)
    if n <= 2 or not n.is_prime():
        raise ValueError("the argument n must be an odd prime")
    Q = quadratic_residues(n)
    Q.remove(0)  # non-zero quad residues
    N = [x for x in srange(1, n) if x not in Q]  # non-zero quad non-residues
    if q not in Q:
        raise ValueError("the order of the finite field must be a quadratic residue modulo n")
    return DuadicCodeOddPair(F, Q, N)

    def iter_indefinite_forms(self) :
        r"""
        Iterate over indices with non-positive discriminant.
        
        TESTS::
        
            sage: from psage.modform.jacobiforms.jacobiformd1nn_fourierexpansion import *
            sage: list(JacobiFormD1NNFilter(2, 2, reduced = False, weak_forms = True).iter_indefinite_forms())
            [(0, -1), (1, 3), (0, 0), (1, -3), (0, 1), (0, -2), (0, 2)]
            sage: list(JacobiFormD1NNFilter(3, 2, reduced = False, weak_forms = True).iter_indefinite_forms())
            [(0, -1), (1, 3), (2, -4), (0, 0), (2, 4), (1, -3), (0, 1), (0, -2), (0, 2)]
            sage: list(JacobiFormD1NNFilter(10, 2, reduced = True, weak_forms = True).iter_indefinite_forms())
            [(0, 0), (0, 1), (0, 2)]
            sage: list(JacobiFormD1NNFilter(10, 3, reduced = True, weak_forms = True).iter_indefinite_forms())
            [(0, 0), (0, 1), (0, 2), (0, 3)]
            sage: list(JacobiFormD1NNFilter(10, 10, reduced = True, weak_forms = True).iter_indefinite_forms())                                                  
            [(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (1, 7), (1, 8), (1, 9), (1, 10), (2, 9), (2, 10)]
        """
        B = self.__bound
        m = self.__m
        fm = Integer(4 * self.__m)
        
        if self.__reduced :
            if self.__weak_forms :
                for n in xrange(0, min(self.__m // 4 + 1, self.__bound)) :
                    for r in xrange( isqrt(fm * n - 1) + 1 if n != 0 else 0, self.__m + 1 ) :
                        yield (n, r)
            else :
                for r in xrange(0, min(self.__m + 1,
                                       isqrt((self.__bound - 1) * fm) + 1) ) :
                    if fm.divides(r**2) :
                        yield (r**2 // fm, r)
        else :
            if self.__weak_forms :
                ## We first determine the reduced indices.
                for n in xrange(0, min(m // 4 + 1, B)) :
                    if n == 0 :
                        r_iteration = range(-m + 1, m + 1)
                    else :
                        r_iteration =   range( -m + 1, -isqrt(fm * n - 1) ) \
                                       + range( isqrt(fm * n - 1) + 1, m + 1 )
                    for r in  r_iteration :
                        for l in range( (- r - isqrt(r**2 - 4 * m * (n - (B - 1))) - 1) // (2 * m) + 1,
                                        (- r + isqrt(r**2 - 4 * m * (n - (B - 1)))) // (2 * m) + 1 ) :
                            if n + l * r + m * l**2 >= B :
                                print l, n, r
                            yield (n + l * r + m * l**2, r + 2 * m * l)
            else :
                if self.__bound > 0 :
                    yield (0,0)

                for n in xrange(1, self.__bound) :
                    if (fm * n).is_square() :
                        rt_fmm = isqrt(fm * n)
                        yield(n, rt_fmm)
                        yield(n, -rt_fmm)
        
        raise StopIteration

    def __init__(self, A, elt=None, check=True):
        """
        TESTS::

            sage: A = FiniteDimensionalAlgebra(GF(3), [Matrix([[1,0], [0,1]]), Matrix([[0,1], [0,0]])])
            sage: A(QQ(4))
            Traceback (most recent call last):
            ...
            TypeError: elt should be a vector, a matrix, or an element of the base field

            sage: B = FiniteDimensionalAlgebra(QQ, [Matrix([[1,0], [0,1]]), Matrix([[0,1], [-1,0]])])
            sage: elt = B(Matrix([[1,1], [-1,1]])); elt
            e0 + e1
            sage: TestSuite(elt).run()
            sage: B(Matrix([[0,1], [1,0]]))
            Traceback (most recent call last):
            ...
            ValueError: matrix does not define an element of the algebra
        """
        AlgebraElement.__init__(self, A)
        k = A.base_ring()
        n = A.degree()
        if elt is None:
            self._vector = vector(k, n)
            self._matrix = Matrix(k, n)
        else:
            if isinstance(elt, int):
                elt = Integer(elt)
            elif isinstance(elt, list):
                elt = vector(elt)
            if A == elt.parent():
                self._vector = elt._vector.base_extend(k)
                self._matrix = elt._matrix.base_extend(k)
            elif k.has_coerce_map_from(elt.parent()):
                e = k(elt)
                if e == 0:
                    self._vector = vector(k, n)
                    self._matrix = Matrix(k, n)
                elif A.is_unitary():
                    self._vector = A._one * e
                    self._matrix = Matrix.identity(k, n) * e
                else:
                    raise TypeError("algebra is not unitary")
            elif is_Vector(elt):
                self._vector = elt.base_extend(k)
                self._matrix = Matrix(k, sum([elt[i] * A.table()[i] for i in xrange(n)]))
            elif is_Matrix(elt):
                if not A.is_unitary():
                    raise TypeError("algebra is not unitary")
                self._vector = A._one * elt
                if not check or sum([self._vector[i]*A.table()[i] for i in xrange(n)]) == elt:
                    self._matrix = elt
                else:
                    raise ValueError("matrix does not define an element of the algebra")
            else:
                raise TypeError("elt should be a vector, a matrix, " +
                                "or an element of the base field")

def _lexi_gen(zeroes=False):
    """
    Return a generator object that produces variable names suitable for the
    generators of a free group.

    INPUT:

    - ``zeroes`` -- Boolean defaulting as ``False``. If ``True``, the
      integers appended to the output string begin at zero at the
      first iteration through the alphabet.

    OUTPUT:

    Python generator object which outputs a character from the alphabet on each
    ``next()`` call in lexicographical order. The integer `i` is appended
    to the output string on the `i^{th}` iteration through the alphabet.

    EXAMPLES::

        sage: from train_track import *
        sage: from train_track.free_group import _lexi_gen
        sage: itr = _lexi_gen()
        sage: F = FreeGroup([next(itr) for i in [1..10]]); F
        Free Group on generators {a, b, c, d, e, f, g, h, i, j}
        sage: it = _lexi_gen()
        sage: [next(it) for i in range(10)]
        ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j']
        sage: itt = _lexi_gen(True)
        sage: [next(itt) for i in range(10)]
        ['a0', 'b0', 'c0', 'd0', 'e0', 'f0', 'g0', 'h0', 'i0', 'j0']
        sage: test = _lexi_gen()
        sage: ls = [next(test) for i in range(3*26)]
        sage: ls[2*26:3*26]
        ['a2', 'b2', 'c2', 'd2', 'e2', 'f2', 'g2', 'h2', 'i2', 'j2', 'k2', 'l2', 'm2',
         'n2', 'o2', 'p2', 'q2', 'r2', 's2', 't2', 'u2', 'v2', 'w2', 'x2', 'y2', 'z2']

    TESTS::

        sage: from train_track import *
        sage: from train_track.free_group import _lexi_gen
        sage: test = _lexi_gen()
        sage: ls = [next(test) for i in range(500)]
        sage: ls[234], ls[260]
        ('a9', 'a10')

    """
    count = Integer(0)
    while True:
        mwrap, ind  = count.quo_rem(26)
        if mwrap == 0 and not(zeroes):
            name = ''
        else:
            name = str(mwrap)
        name = chr(ord('a') + ind) + name
        yield name
        count = count + 1

    def __init__(self, q, names="a",  modulus=None, repr="poly"):
        """
        Initialize ``self``.

        TESTS::

            sage: k.<a> = GF(2^100, modulus='strangeinput')
            Traceback (most recent call last):
            ...
            ValueError: Modulus parameter not understood
            sage: k.<a> = GF(2^20) ; type(k)
            <class 'sage.rings.finite_rings.finite_field_ntl_gf2e.FiniteField_ntl_gf2e_with_category'>
            sage: loads(dumps(k)) is k
            True
            sage: k1.<a> = GF(2^16)
            sage: k2.<a> = GF(2^17)
            sage: k1 == k2
            False
            sage: k3 = k1._finite_field_ext_pari_()
            sage: k1 == k3
            False

            sage: TestSuite(k).run()

            sage: k.<a> = GF(2^64)
            sage: k._repr_option('element_is_atomic')
            False
            sage: P.<x> = PolynomialRing(k)
            sage: (a+1)*x # indirect doctest
            (a + 1)*x
        """
        late_import()
        q = Integer(q)
        if q < 2:
            raise ValueError("q must be a 2-power")
        k = q.exact_log(2)
        if q != 1 << k:
            raise ValueError("q must be a 2-power")
        p = Integer(2)
        FiniteField.__init__(self, GF(p), names, normalize=True)

        self._kwargs = {'repr':repr}
        self._is_conway = False

        if modulus is None or modulus == 'default':
            if exists_conway_polynomial(p, k):
                modulus = "conway"
            else:
                modulus = "minimal_weight"
        if modulus == "conway":
            modulus = conway_polynomial(p, k)
            self._is_conway = True
        if is_Polynomial(modulus):
            modulus = modulus.list()
        self._cache = Cache_ntl_gf2e(self, k, modulus)
        self._polynomial = {}

    def __init__(self, q, name="a", modulus=None, repr="poly", cache=False):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: k.<a> = GF(2^3)
            sage: j.<b> = GF(3^4)
            sage: k == j
            False

            sage: GF(2^3,'a') == copy(GF(2^3,'a'))
            True
            sage: TestSuite(GF(2^3, 'a')).run()
        """
        self._kwargs = {}

        if repr not in ["int", "log", "poly"]:
            raise ValueError("Unknown representation %s" % repr)

        q = Integer(q)
        if q < 2:
            raise ValueError("q must be a prime power")
        F = q.factor()
        if len(F) > 1:
            raise ValueError("q must be a prime power")
        p = F[0][0]
        k = F[0][1]

        if q >= 1 << 16:
            raise ValueError("q must be < 2^16")

        from .finite_field_constructor import GF

        FiniteField.__init__(self, GF(p), name, normalize=False)

        self._kwargs["repr"] = repr
        self._kwargs["cache"] = cache

        from sage.rings.polynomial.polynomial_element import is_Polynomial

        if not is_Polynomial(modulus):
            from sage.misc.superseded import deprecation

            deprecation(
                16930,
                "constructing a FiniteField_givaro without giving a polynomial as modulus is deprecated, use the more general FiniteField constructor instead",
            )
            R = GF(p)["x"]
            if modulus is None or isinstance(modulus, str):
                modulus = R.irreducible_element(k, algorithm=modulus)
            else:
                modulus = R(modulus)

        self._cache = Cache_givaro(self, p, k, modulus, repr, cache)
        self._modulus = modulus

    def __classcall_private__(cls, begin, end=None, step=Integer(1), middle_point=None):
        """
        TESTS::

            sage: IntegerRange(2,5,0)
            Traceback (most recent call last):
            ...
            ValueError: IntegerRange() step argument must not be zero
            sage: IntegerRange(2) is IntegerRange(0, 2)
            True
            sage: IntegerRange(1.0)
            Traceback (most recent call last):
            ...
            TypeError: end must be Integer or Infinity, not <type 'sage.rings.real_mpfr.RealLiteral'>
        """
        if isinstance(begin, int): begin = Integer(begin)
        if isinstance(end, int): end = Integer(end)
        if isinstance(step,int): step = Integer(step)

        if end is None:
            end = begin
            begin = Integer(0)
        # check of the arguments
        if not isinstance(begin, (Integer, MinusInfinity, PlusInfinity)):
            raise TypeError("begin must be Integer or Infinity, not %r" % type(begin))
        if not isinstance(end, (Integer, MinusInfinity, PlusInfinity)):
            raise TypeError("end must be Integer or Infinity, not %r" % type(end))
        if not isinstance(step, Integer):
            raise TypeError("step must be Integer, not %r" % type(step))
        if step.is_zero():
            raise ValueError("IntegerRange() step argument must not be zero")

        # If begin and end are infinite, middle_point and step will defined the set.
        if begin == -Infinity and end == Infinity:
            if middle_point is None:
                raise ValueError("Can't iterate over this set, please provide middle_point")

        # If we have a middle point, we go on the special enumeration way...
        if middle_point is not None:
            return IntegerRangeFromMiddle(begin, end, step, middle_point)

        if (begin == -Infinity) or (begin == Infinity):
            raise ValueError("Can't iterate over this set: It is impossible to begin an enumeration with plus/minus Infinity")

        # Check for empty sets
        if step > 0 and begin >= end or step < 0 and begin <= end:
            return IntegerRangeEmpty()

        if end != Infinity and end != -Infinity:
            # Normalize the input
            sgn = 1 if step > 0 else -1
            end = begin+((end-begin-sgn)//(step)+1)*step
            return IntegerRangeFinite(begin, end, step)
        else:
            return IntegerRangeInfinite(begin, step)

    def __init__(self, q, names="a", modulus=None, repr="poly"):
        """
        Initialize ``self``.

        TESTS::

            sage: k.<a> = GF(2^100, modulus='strangeinput')
            Traceback (most recent call last):
            ...
            ValueError: no such algorithm for finding an irreducible polynomial: strangeinput
            sage: k.<a> = GF(2^20) ; type(k)
            <class 'sage.rings.finite_rings.finite_field_ntl_gf2e.FiniteField_ntl_gf2e_with_category'>
            sage: loads(dumps(k)) is k
            True
            sage: k1.<a> = GF(2^16)
            sage: k2.<a> = GF(2^17)
            sage: k1 == k2
            False
            sage: k3.<a> = GF(2^16, impl="pari_ffelt")
            sage: k1 == k3
            False

            sage: TestSuite(k).run()

            sage: k.<a> = GF(2^64)
            sage: k._repr_option('element_is_atomic')
            False
            sage: P.<x> = PolynomialRing(k)
            sage: (a+1)*x # indirect doctest
            (a + 1)*x
        """
        late_import()
        q = Integer(q)
        if q < 2:
            raise ValueError("q must be a 2-power")
        k = q.exact_log(2)
        if q != 1 << k:
            raise ValueError("q must be a 2-power")
        FiniteField.__init__(self, GF2, names, normalize=True)

        self._kwargs = {'repr':repr}

        from sage.rings.polynomial.polynomial_element import is_Polynomial
        if not is_Polynomial(modulus):
            from sage.misc.superseded import deprecation
            deprecation(16983, "constructing a FiniteField_ntl_gf2e without giving a polynomial as modulus is deprecated, use the more general FiniteField constructor instead")
            R = GF2['x']
            if modulus is None or isinstance(modulus, str):
                modulus = R.irreducible_element(k, algorithm=modulus)
            else:
                modulus = R(modulus)

        self._cache = Cache_ntl_gf2e(self, k, modulus)
        self._modulus = modulus

    def __init__(self, q, name="a",  modulus=None, repr="poly", cache=False):
        """
        Initialize ``self``.

        EXAMPLES::

            sage: k.<a> = GF(2^3)
            sage: j.<b> = GF(3^4)
            sage: k == j
            False

            sage: GF(2^3,'a') == copy(GF(2^3,'a'))
            True
            sage: TestSuite(GF(2^3, 'a')).run()
        """
        self._kwargs = {}

        if repr not in ['int', 'log', 'poly']:
            raise ValueError("Unknown representation %s"%repr)

        q = Integer(q)
        if q < 2:
            raise ValueError("q  must be a prime power")
        F = q.factor()
        if len(F) > 1:
            raise ValueError("q must be a prime power")
        p = F[0][0]
        k = F[0][1]

        if q >= 1<<16:
            raise ValueError("q must be < 2^16")

        import constructor
        FiniteField.__init__(self, constructor.FiniteField(p), name, normalize=False)

        self._kwargs['repr'] = repr
        self._kwargs['cache'] = cache

        if modulus is None or modulus == 'conway':
            if k == 1:
                modulus = 'random' # this will use the gfq_factory_pk function.
            elif ConwayPolynomials().has_polynomial(p, k):
                from sage.rings.finite_rings.conway_polynomials import conway_polynomial
                modulus = conway_polynomial(p, k)
            elif modulus is None:
                modulus = 'random'
            else:
                raise ValueError("Conway polynomial not found")

        from sage.rings.polynomial.polynomial_element import is_Polynomial
        if is_Polynomial(modulus):
            modulus = modulus.list()

        self._cache = Cache_givaro(self, p, k, modulus, repr, cache)

def _big_delta_coeff(aa, bb, cc, prec=None):
    r"""
    Calculates the Delta coefficient of the 3 angular momenta for
    Racah symbols. Also checks that the differences are of integer
    value.

    INPUT:

    -  ``aa`` - first angular momentum, integer or half integer

    -  ``bb`` - second angular momentum, integer or half integer

    -  ``cc`` - third angular momentum, integer or half integer

    -  ``prec`` - precision of the ``sqrt()`` calculation

    OUTPUT:

    double - Value of the Delta coefficient

    EXAMPLES::

        sage: from sage.functions.wigner import _big_delta_coeff
        sage: _big_delta_coeff(1,1,1)
        1/2*sqrt(1/6)
    """
    if int(aa + bb - cc) != (aa + bb - cc):
        raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
    if int(aa + cc - bb) != (aa + cc - bb):
        raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
    if int(bb + cc - aa) != (bb + cc - aa):
        raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
    if (aa + bb - cc) < 0:
        return 0
    if (aa + cc - bb) < 0:
        return 0
    if (bb + cc - aa) < 0:
        return 0

    maxfact = max(aa + bb - cc, aa + cc - bb, bb + cc - aa, aa + bb + cc + 1)
    _calc_factlist(maxfact)

    argsqrt = Integer(_Factlist[int(aa + bb - cc)] * \
                          _Factlist[int(aa + cc - bb)] * \
                          _Factlist[int(bb + cc - aa)]) / \
                          Integer(_Factlist[int(aa + bb + cc + 1)])

    ressqrt = argsqrt.sqrt(prec)
    if isinstance(ressqrt, ComplexNumber):
        res = ressqrt.real()
    else:
        res = ressqrt
    return res

def _ab( m):
    r"""
    OUTPUT
        The unique pair of integers $a,b$ such that $m=a^2b$
        with squarefree $b$.

    INPUT
        m -- a nonzero integer
    """
    m = Integer(m)
    b = m.squarefree_part()
    a = (m//b).isqrt()
    return a, b

    def __getitem__(self, level):
        """
        Return the modular polynomial of given level, or an error if
        there is no such polynomial in the database.

        EXAMPLES::

            sage: DBMP = ClassicalModularPolynomialDatabase()
            sage: f = DBMP[29]                                 # optional - database_kohel
            sage: f.degree()                                   # optional - database_kohel
            58
            sage: f.coefficient([28,28])                       # optional - database_kohel
            400152899204646997840260839128

            sage: DBMP[50]                                     # optional - database_kohel
            Traceback (most recent call last):
            ...
            LookupError: filename .../kohel/PolMod/Cls/pol.050.dbz does not exist
        """
        from sage.rings.integer import Integer
        from sage.rings.integer_ring import IntegerRing
        from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

        if self.model in ("Atk","Eta"):
            level = Integer(level)
            if not level.is_prime():
                raise TypeError("Argument level (= %s) must be prime."%level)
        elif self.model in ("AtkCrr","EtaCrr"):
            N = Integer(level[0])
            if not N in (2,3,5,7,13):
                raise TypeError("Argument level (= %s) must be prime."%N)
        modpol = self._dbpath(level)
        coeff_list = _dbz_to_integer_list(modpol)
        if self.model == "Cls":
            P = PolynomialRing(IntegerRing(),2,"j")
        else:
            P = PolynomialRing(IntegerRing(),2,"x,j")
        poly = {}
        if self.model == "Cls":
            if level == 1:
                return P({(1,0):1,(0,1):-1})
            for cff in coeff_list:
                i = cff[0]
                j = cff[1]
                poly[(i,j)] = Integer(cff[2])
                if i != j:
                    poly[(j,i)] = Integer(cff[2])
        else:
            for cff in coeff_list:
                poly[(cff[0],cff[1])] = Integer(cff[2])
        return P(poly)

    def __getitem__(self,level):
        """
        Return the modular polynomial of given level, or an error if
        there is no such polynomial in the database.

        EXAMPLES:
            sage: DBMP = ClassicalModularPolynomialDatabase()  # optional - database_kohel
            sage: f = DBMP[29]                                 #optional
            sage: f.degree()                                   #optional
            58
            sage: f.coefficient([28,28])                       #optional
            400152899204646997840260839128

            sage: DBMP[50]                                     #optional
            Traceback (most recent call last):
            ...
            RuntimeError: No database entry for modular polynomial of level 50
        """
        if self.model in ("Atk","Eta"):
            level = Integer(level)
            if not level.is_prime():
                raise TypeError("Argument level (= %s) must be prime."%level)
        elif self.model in ("AtkCrr","EtaCrr"):
            N = Integer(level[0])
            if not N in (2,3,5,7,13):
                raise TypeError("Argument level (= %s) must be prime."%N)
        modpol = self._dbpath(level)
        try:
            coeff_list = _dbz_to_integer_list(modpol)
        except RuntimeError as msg:
            print(msg)
            raise RuntimeError("No database entry for modular polynomial of level %s"%level)
        if self.model == "Cls":
            P = PolynomialRing(IntegerRing(),2,"j")
        else:
            P = PolynomialRing(IntegerRing(),2,"x,j")
        poly = {}
        if self.model == "Cls":
            if level == 1:
                return P({(1,0):1,(0,1):-1})
            for cff in coeff_list:
                i = cff[0]
                j = cff[1]
                poly[(i,j)] = Integer(cff[2])
                if i != j:
                    poly[(j,i)] = Integer(cff[2])
        else:
            for cff in coeff_list:
                poly[(cff[0],cff[1])] = Integer(cff[2])
        return P(poly)