def Pall_mass_density_at_odd_prime(self, p):
    """
    Returns the local representation density of a form (for
    representing itself) defined over `ZZ`, at some prime `p>2`.

    REFERENCES:
        Pall's article "The Weight of a Genus of Positive n-ary Quadratic Forms"
        appearing in Proc. Symp. Pure Math. VIII (1965), pp95--105.

    INPUT:
        `p` -- a prime number > 2.

    OUTPUT:
        a rational number.

    EXAMPLES::

        sage: Q = QuadraticForm(ZZ, 3, [1,0,0,1,0,1])
        sage: Q.Pall_mass_density_at_odd_prime(3)
        [(0, Quadratic form in 3 variables over Integer Ring with coefficients:
        [ 1 0 0 ]
        [ * 1 0 ]
        [ * * 1 ])] [(0, 3, 8)] [8/9] 8/9
        8/9
    """
    ## Check that p is a positive prime -- unnecessary since it's done implicitly in the next step. =)
    if p<=2:
        raise TypeError, "Oops!  We need p to be a prime > 2."

    ## Step 1: Obtain a p-adic (diagonal) local normal form, and
    ## compute the invariants for each Jordan block.
    jordan_list = self.jordan_blocks_by_scale_and_unimodular(p)
    modified_jordan_list = [(a, Q.dim(), Q.det())  for (a,Q) in jordan_list]     ## List of pairs (scale, det)
    #print jordan_list
    #print modified_jordan_list

    ## Step 2: Compute the list of local masses for each Jordan block
    jordan_mass_list = []
    for (s,n,d) in modified_jordan_list:
        generic_factor = prod([1 - p**(-2*j)  for j in range(1, floor((n-1)/2)+1)])
        #print "generic factor: ", generic_factor
        if (n % 2 == 0):
            m = n/2
            generic_factor *= (1 + legendre_symbol(((-1)**m) * d, p) * p**(-m))
        #print "jordan_mass: ", generic_factor
        jordan_mass_list = jordan_mass_list + [generic_factor]

    ## Step 3: Compute the local mass $\al_p$ at p.
        MJL = modified_jordan_list
    s = len(modified_jordan_list)
    M = [sum([MJL[j][1]  for j in range(i, s)])  for i in range(s-1)]    ## Note: It's s-1 since we don't need the last M.
    #print "M = ", M
    nu = sum([M[i] * MJL[i][0] * MJL[i][1]  for i in range(s-1)]) - ZZ(sum([J[0] * J[1] * (J[1]-1)  for J in MJL]))/ZZ(2)
    p_mass = prod(jordan_mass_list)
    p_mass *= 2**(s-1) * p**nu

    print jordan_list, MJL, jordan_mass_list, p_mass

    ## Return the result
    return p_mass

    def number_of_classes(self, invertible=False, q=None):
        """
        Return the number of similarity classes of matrices of type ``self``.

        IMPUT:

        - ``invertible`` -- Boolean; return number of invertible classes if set
          to ``True``

        - ``q`` -- An integer or an indeterminate

        EXAMPLES::

            sage: tau = SimilarityClassType([[1, [1]], [1, [1]]])
            sage: tau.number_of_classes()
            1/2*q^2 - 1/2*q
        """
        if q is None:
            q = ZZ["q"].gen()
        if self.size() == 0:
            return q.parent().one()
        list_of_degrees = [PT.degree() for PT in self]
        maximum_degree = max(list_of_degrees)
        numerator = prod(
            [
                prod([primitives(d + 1, invertible=invertible, q=q) - i for i in range(list_of_degrees.count(d + 1))])
                for d in range(maximum_degree)
            ]
        )
        tau_list = list(self)
        D = dict((i, tau_list.count(i)) for i in tau_list)
        denominator = reduce(mul, [factorial(D[primary_type]) for primary_type in D.keys()])
        return numerator / denominator

def q_binomial(n,k,p=None):
    """
    Returns the ``q``-binomial coefficient.

    If ``p`` is unspecified, then it defaults to using the generator ``q`` for
    a univariate polynomial ring over the integers.

    EXAMPLES::

        sage: import sage.combinat.q_analogues as q_analogues
        sage: q_analogues.q_binomial(4,2)
        q^4 + q^3 + 2*q^2 + q + 1
        sage: p = ZZ['p'].0
        sage: q_analogues.q_binomial(4,2,p)
        p^4 + p^3 + 2*p^2 + p + 1

    The ``q``-analogue of ``binomial(n,k)`` is currently only defined for
    ``n`` a nonnegative integer, it is zero for negative k  (trac #11411)::

        sage: q_analogues.q_binomial(5, -1)
        0
    """
    if not (n in ZZ and k in ZZ):
        raise ValueError, "Argument (%s, %s) must be integers."%(n, k)
    if n < 0:
        raise NotImplementedError
    if 0 <= k and k <= n:
        k=min(k, n-k)
        return (prod(q_int(j, p) for j in range(n-k+1, n+1)) /
                prod(q_int(j, p) for j in range(1, k+1)))
    else:
        return 0

def q_catalan_number(n,p=None):
    """
    Returns the `q`-Catalan number of index `n`.

    If `p` is unspecified, then it defaults to using the generator `q` for
    a univariate polynomial ring over the integers.

    There are several `q`-Catalan numbers. This procedure
    returns the one which can be written using the `q`-binomial coefficients.

    EXAMPLES::

        sage: from sage.combinat.q_analogues import q_catalan_number
        sage: q_catalan_number(4)
        q^12 + q^10 + q^9 + 2*q^8 + q^7 + 2*q^6 + q^5 + 2*q^4 + q^3 + q^2 + 1
        sage: p = ZZ['p'].0
        sage: q_catalan_number(4,p)
        p^12 + p^10 + p^9 + 2*p^8 + p^7 + 2*p^6 + p^5 + 2*p^4 + p^3 + p^2 + 1

    The `q`-Catalan number of index `n` is only defined for `n` a
    nonnegative integer (:trac:`11411`)::

        sage: q_catalan_number(-2)
        Traceback (most recent call last):
        ...
        ValueError: Argument (-2) must be a nonnegative integer.
    """
    if n in ZZ and n >= 0:
        return prod(q_int(j, p) for j in range(n+2, 2*n+1)) / prod(q_int(j, p) for j in range(2,n+1))
    else:
        raise ValueError("Argument (%s) must be a nonnegative integer." %n)

def DuadicCodeOddPair(F,S1,S2):
    """
    Constructs the "odd pair" of duadic codes associated to the
    "splitting" S1, S2 of n.

    .. warning::

       Maybe the splitting should be associated to a sum of
       q-cyclotomic cosets mod n, where q is a *prime*.

    EXAMPLES::

        sage: from sage.coding.code_constructions import is_a_splitting
        sage: n = 11; q = 3
        sage: C = cyclotomic_cosets(q,n); C
        [[0], [1, 3, 4, 5, 9], [2, 6, 7, 8, 10]]
        sage: S1 = C[1]
        sage: S2 = C[2]
        sage: is_a_splitting(S1,S2,11)
        (True, 2)
        sage: DuadicCodeOddPair(GF(q),S1,S2)
        (Linear code of length 11, dimension 6 over Finite Field of size 3,
         Linear code of length 11, dimension 6 over Finite Field of size 3)

    This is consistent with Theorem 6.1.3 in [HP]_.
    """
    n = max(S1+S2)+1
    if not(is_a_splitting(S1,S2,n)):
        raise TypeError, "%s, %s must be a splitting of %s."%(S1,S2,n)
    q = F.order()
    k = Mod(q,n).multiplicative_order()
    FF = GF(q**k,"z")
    z = FF.gen()
    zeta = z**((q**k-1)/n)
    P1 = PolynomialRing(FF,"x")
    x = P1.gen()
    g1 = prod([x-zeta**i for i in S1+[0]])
    g2 = prod([x-zeta**i for i in S2+[0]])
    j = sum([x**i/n for i in range(n)])
    P2 = PolynomialRing(F,"x")
    x = P2.gen()
    coeffs1 = [lift2smallest_field(c)[0] for c in (g1+j).coeffs()]
    coeffs2 = [lift2smallest_field(c)[0] for c in (g2+j).coeffs()]
    gg1 = P2(coeffs1)
    gg2 = P2(coeffs2)
    C1 = CyclicCodeFromGeneratingPolynomial(n,gg1)
    C2 = CyclicCodeFromGeneratingPolynomial(n,gg2)
    return C1,C2

 def _gap_init_(self):
     r"""
     Return string that defines corresponding abelian group in GAP.
     
     EXAMPLES::
     
         sage: G = AbelianGroup([2,3,9])
         sage: G._gap_init_()
         'AbelianGroup([2, 3, 9])'
         sage: gap(G)
         Group( [ f1, f2, f3 ] )
     
     Only works for finite groups.
     
     ::
     
         sage: G = AbelianGroup(3,[0,3,4],names="abc"); G
         Multiplicative Abelian Group isomorphic to Z x C3 x C4
         sage: G._gap_init_()
         Traceback (most recent call last):
         ...
         TypeError: abelian groups in GAP are finite, but self is infinite
     """
     # TODO: Use the package polycyclic has AbelianPcpGroup, which can handle
     # the infinite case but it is a GAP package not GPL'd.
     # Use this when the group is infinite...
     if (False and prod(self.invariants())==0):   # if False for now...
         return 'AbelianPcpGroup(%s)'%self.invariants()
     if not self.is_finite():
         raise TypeError, "abelian groups in GAP are finite, but self is infinite"
     return 'AbelianGroup(%s)'%self.invariants()

    def nu3(self):
        r"""
        Return the number of elliptic points of order 3 for this congruence
        subgroup `\Gamma_0(N)`. The number of these is given by a standard formula:
        0 if `N` is divisible by 9 or any prime congruent to -1 mod 3, and
        otherwise `2^d` where d is the number of primes other than 3 dividing `N`.

        EXAMPLE::

            sage: Gamma0(2).nu3()
            0
            sage: Gamma0(3).nu3()
            1
            sage: Gamma0(9).nu3()
            0
            sage: Gamma0(7).nu3()
            2
            sage: Gamma0(21).nu3()
            2
            sage: Gamma0(1729).nu3()
            8
        """
        n = self.level()
        if (n % 9 == 0):
            return ZZ(0)
        return prod([ 1 + kronecker_symbol(-3, p) for p, _ in n.factor()])

    def cardinality(self):
        r"""
        Return the cardinality of ``self``.

        The number of ordered set partitions of a set of length `k` with
        composition shape `\mu` is equal to

        .. MATH::

            \frac{k!}{\prod_{\mu_i \neq 0} \mu_i!}.

        EXAMPLES::

            sage: OrderedSetPartitions(5,[2,3]).cardinality()
            10
            sage: OrderedSetPartitions(0, []).cardinality()
            1
            sage: OrderedSetPartitions(0, [0]).cardinality()
            1
            sage: OrderedSetPartitions(0, [0,0]).cardinality()
            1
            sage: OrderedSetPartitions(5, [2,0,3]).cardinality()
            10
        """
        return factorial(len(self._set))/prod([factorial(i) for i in self.c])

    def statistic(self, func, q=None):
        """
        Return

        .. MATH::

            prod_{(d, \lambda)\in \tau} n_\lambda(q^d)

        where `n_\lambda(q)` is the value returned by ``func`` on the input
        `\lambda`.

        INPUT:

        - ``func`` -- a function that takes a partition to a polynomial in ``q``

        - ``q`` -- an integer or an indeterminate

        EXAMPLES::

            sage: tau = SimilarityClassType([[1, [1]], [1, [2, 1]], [2, [1, 1]]])
            sage: from sage.combinat.similarity_class_type import fq
            sage: tau.statistic(lambda la: prod([fq(m) for m in la.to_exp()]))
            (q^9 - 3*q^8 + 2*q^7 + 2*q^6 - 4*q^5 + 4*q^4 - 2*q^3 - 2*q^2 + 3*q - 1)/q^9
            sage: q = ZZ['q'].gen()
            sage: tau.statistic(lambda la: q**la.size(), q = q)
            q^8
        """
        if q is None:
            q = FractionField(ZZ["q"]).gen()
        return prod([PT.statistic(func, q=q) for PT in self])

def centralizer_group_cardinality(la, q=None):
    r"""
    Return the cardinality of the centralizer group in `GL_n(\mathbf{F}_q)`
    of a nilpotent matrix whose Jordan blocks are given by ``la``.

    INPUT:

    - ``lambda`` -- a partition

    - ``q`` -- an integer or an indeterminate

    OUTPUT:

    A polynomial function of ``q``.

    EXAMPLES::

        sage: from sage.combinat.similarity_class_type import centralizer_group_cardinality
        sage: q = ZZ['q'].gen()
        sage: centralizer_group_cardinality(Partition([2, 1]))
        q^5 - 2*q^4 + q^3
    """
    if q is None:
        q = ZZ["q"].gen()
    return q ** centralizer_algebra_dim(la) * prod([fq(m, q=q) for m in la.to_exp()])

    def _Chow_group_free(self):
        r"""
        Return the relations coming from the free part of the Chow group

        OUTPUT:

        A tuple containing the elements of $Hom(A_{d-1,\text{free}},
        F^\times)$, including the identity.

        EXAMPLES::

            sage: fan = NormalFan(ReflexivePolytope(2, 0))
            sage: X = ToricVariety(fan, base_field=GF(7))
            sage: X.Chow_group().degree(1)
            C3 x Z
            sage: enum = X.point_set()._naive_enumerator()
            sage: enum._Chow_group_free()
            ((1, 1, 1), (2, 2, 2), (3, 3, 3), (4, 4, 4), (5, 5, 5), (6, 6, 6))
        """
        units = self.units()
        result = []
        rays = self.fan().rays() + self.fan().virtual_rays()
        ker = rays.matrix().integer_kernel().matrix()
        for phases in CartesianProduct(*([units] * ker.nrows())):
            phases = tuple(prod(mu**exponent for mu, exponent in zip(phases, column))
                           for column in ker.columns())
            result.append(phases)
        return tuple(sorted(result))

    def exp(self, exponents):
        r"""
        Return unit with given exponents with respect to group generators.

        INPUT:

        - ``u`` -- Any object from which an element of the unit
          group's number field `K` may be constructed; an error is
          raised if an element of `K` cannot be constructed from u, or
          if the element constructed is not a unit.

        OUTPUT: a list of integers giving the exponents of ``u`` with
        respect to the unit group's basis.
        
        EXAMPLES::

            sage: x = polygen(QQ)
            sage: K.<z> = CyclotomicField(13)
            sage: UK = UnitGroup(K)
            sage: [UK.log(u) for u in UK.gens()]
            [(1, 0, 0, 0, 0, 0),
             (0, 1, 0, 0, 0, 0),
             (0, 0, 1, 0, 0, 0),
             (0, 0, 0, 1, 0, 0),
             (0, 0, 0, 0, 1, 0),
             (0, 0, 0, 0, 0, 1)]
            sage: vec = [65,6,7,8,9,10]
            sage: unit = UK.exp(vec)
            sage: UK.log(unit)
            (13, 6, 7, 8, 9, 10)
            sage: UK.exp(UK.log(u)) == u.value()
            True
        """
        return prod([u**e for u,e in zip(self.gens_values(),exponents)], self.number_field().one_element())

 def __call__( self, g ):
     """
     Some python code for wrapping GAP's Images function but only for
     permutation groups. Returns an error if g is not in G.
     
     EXAMPLES::
     
         sage: H = AbelianGroup(3, [2,3,4], names="abc")
         sage: a,b,c = H.gens()
         sage: G = AbelianGroup(2, [2,3], names="xy")
         sage: x,y = G.gens()
         sage: phi = AbelianGroupMorphism(G,H,[x,y],[a,b])
         sage: phi(y*x)
         a*b
         sage: phi(y^2)
         b^2
     """
     G = g.parent()
     w = g.word_problem(self.domaingens)
     n = len(w)
     #print w,g.word_problem(self.domaingens)
     # g.word_problem is faster in general than word_problem(g)
     gens = self.codomaingens
     h = prod([gens[(self.domaingens).index(w[i][0])]**(w[i][1]) for i in range(n)])
     return h

def Ht(mu, q=None, t=None, pi=None):
    """
    Returns the symmetric Macdonald polynomial using the Haiman,
    Haglund, and Loehr formula.

    Note that if both `q` and `t` are specified, then they must have the
    same parent.

    REFERENCE:

    - J. Haglund, M. Haiman, N. Loehr.
      *A combinatorial formula for non-symmetric Macdonald polynomials*.
      :arXiv:`math/0601693v3`.

    EXAMPLES::

        sage: from sage.combinat.sf.ns_macdonald import Ht
        sage: HHt = SymmetricFunctions(QQ['q','t'].fraction_field()).macdonald().Ht()
        sage: Ht([0,0,1])
        x0 + x1 + x2
        sage: HHt([1]).expand(3)
        x0 + x1 + x2
        sage: Ht([0,0,2])
        x0^2 + (q + 1)*x0*x1 + x1^2 + (q + 1)*x0*x2 + (q + 1)*x1*x2 + x2^2
        sage: HHt([2]).expand(3)
        x0^2 + (q + 1)*x0*x1 + x1^2 + (q + 1)*x0*x2 + (q + 1)*x1*x2 + x2^2
    """
    P, q, t, n, R, x = _check_muqt(mu, q, t, pi)
    res = 0
    for a in n:
        weight = a.weight()
        res += q**a.maj()*t**a.inv()*prod( x[i]**weight[i] for i in range(len(weight)) )
    return res

 def __pow__(self, n):
     """
     requires that len(invs) = n
     """
     if not isinstance(n, (int, long, Integer)):
         raise TypeError, "Argument n (= %s) must be an integer."%n
     n = int(n)
     M = self.parent()
     N = M.ngens()
     invs = M.invariants()
     if n < 0:  
         L =[n*self.list()[i]%M.gen(i).order() for i in range(M.ngens())]
         return prod([M.gen(i)**L[i] for i in range(M.ngens())])
         #m = LCM(invs) ## Not very efficient version
         #pw = (n)%m
         #x = self**pw
         #return x
     elif n == 0:
         return M(1)
     elif n == 1:
         return self
     elif n == 2:
         return self * self
     k = n//2
     return self**k * self**(n-k)

def _monomial_exponent_to_lower_factorial(me, x):
    r"""
    Converts a tuple of exponents to the monomial obtained by replacing
    each me[i] with `x_i*(x_i - 1)*\cdots*(x_i - a_i + 1)`
    
    EXAMPLES::
    
        sage: from sage.combinat.misc import _monomial_exponent_to_lower_factorial
        sage: R.<x,y,z> = QQ[]
        sage: a = R.gens()
        sage: _monomial_exponent_to_lower_factorial(([1,0,0]),a)
        x
        sage: _monomial_exponent_to_lower_factorial(([2,0,0]),a)
        x^2 - x
        sage: _monomial_exponent_to_lower_factorial(([0,2,0]),a)
        y^2 - y
        sage: _monomial_exponent_to_lower_factorial(([1,1,0]),a)
        x*y
        sage: _monomial_exponent_to_lower_factorial(([1,1,2]),a)
        x*y*z^2 - x*y*z
        sage: _monomial_exponent_to_lower_factorial(([2,2,2]),a)
        x^2*y^2*z^2 - x^2*y^2*z - x^2*y*z^2 - x*y^2*z^2 + x^2*y*z + x*y^2*z + x*y*z^2 - x*y*z
    """
    terms = []
    for i in range(len(me)):
        for j in range(me[i]):
            terms.append(x[i] - j)
    return prod(terms)

    def __getitem__(self, support):
        r"""
        INPUT:

         - ``support`` - a proper subset of the index_set, as a list or set.

        Returns the cyclicaly decreasing element associated with ``support``.

        EXAMPLES::
:
            sage: W = WeylGroup(["A", 5, 1])
            sage: W.pieri_factors()[[0,1,2,3,5]].reduced_word()
            [3, 2, 1, 0, 5]
            sage: W.pieri_factors()[[0,1,3,4,5]].reduced_word()
            [1, 0, 5, 4, 3]
            sage: W.pieri_factors()[[0,1,2,3,4]].reduced_word()
            [4, 3, 2, 1, 0]

        """
        index_set = sorted(self.W.index_set())
        support   = sorted(support)
        assert set(support).issubset(set(index_set))
        assert support != index_set
        if len(support) == 0:
            return self.W.one()
        s = self.W.simple_reflections()
        i = 0
        while i < len(support) and support[i] == index_set[i]:
            i += 1
        # This finds the first hole: either ley[i] is maximal or support[i] < support[i+1]+1
        return prod((s[j] for j in list(reversed(support[0:i])) + list(reversed(support[i:]))), self.W.one())

    def cardinality(self):
        """
        Returns the number of integer necklaces with the evaluation e.
        
        EXAMPLES::
        
            sage: Necklaces([]).cardinality()
            0
            sage: Necklaces([2,2]).cardinality()
            2
            sage: Necklaces([2,3,2]).cardinality()
            30
        
        Check to make sure that the count matches up with the number of
        Lyndon words generated.
        
        ::
        
            sage: comps = [[],[2,2],[3,2,7],[4,2]]+Compositions(4).list()
            sage: ns = [ Necklaces(comp) for comp in comps]
            sage: all( [ n.cardinality() == len(n.list()) for n in ns] )
            True
        """
        evaluation = self.e
        le = list(evaluation)
        if len(le) == 0:
            return 0

        n = sum(le)

        return sum([euler_phi(j)*factorial(n/j) / prod([factorial(ni/j) for ni in evaluation]) for j in divisors(gcd(le))])/n

    def cardinality(self):
        """
        Returns the number of Lyndon words with the evaluation e.

        EXAMPLES::

            sage: LyndonWords([]).cardinality()
            0
            sage: LyndonWords([2,2]).cardinality()
            1
            sage: LyndonWords([2,3,2]).cardinality()
            30

        Check to make sure that the count matches up with the number of
        Lyndon words generated.

        ::

            sage: comps = [[],[2,2],[3,2,7],[4,2]]+Compositions(4).list()
            sage: lws = [ LyndonWords(comp) for comp in comps]
            sage: all( [ lw.cardinality() == len(lw.list()) for lw in lws] )
            True
        """
        evaluation = self.e
        le = __builtin__.list(evaluation)
        if len(evaluation) == 0:
            return 0

        n = sum(evaluation)

        return sum([moebius(j)*factorial(n/j) / prod([factorial(ni/j) for ni in evaluation]) for j in divisors(gcd(le))])/n

def psi(N):
    """
    The index `[\Gamma : \Gamma_0(N)]`, where `\Gamma = GL(2, R)` for `R` the
    corresponding ring of integers, and `\Gamma_0(N)` standard congruence
    subgroup.

    EXAMPLES::

        sage: from sage.modular.modsym.p1list_nf import psi
        sage: k.<a> = NumberField(x^2 + 23)
        sage: N = k.ideal(3, a - 1)
        sage: psi(N)
        4

    ::

        sage: k.<a> = NumberField(x^2 + 23)
        sage: N = k.ideal(5)
        sage: psi(N)
        26
    """
    if not N.is_integral():
        raise ValueError("psi only defined for integral ideals")

    from sage.misc.misc import prod
    return prod([(np+1)*np**(e-1) \
                     for np,e in [(p.absolute_norm(),e) \
                                  for p,e in N.factor()]])