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

    If `q` 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, q) for j in range(n+2, 2*n+1)) / prod(q_int(j, q) for j in range(2,n+1))
    else:
        raise ValueError("Argument (%s) must be a nonnegative integer." %n)

        def weyl_dimension(self, highest_weight):
            r"""
            Return the dimension of the highest weight representation of highest weight ``highest_weight``.

            EXAMPLES::

                sage: RootSystem(['A',3]).ambient_lattice().weyl_dimension([2,1,0,0])
                20
                sage: P = RootSystem(['C',2]).weight_lattice()
                sage: La = P.basis()
                sage: P.weyl_dimension(La[1]+La[2])
                16

                sage: type(RootSystem(['A',3]).ambient_lattice().weyl_dimension([2,1,0,0]))
                <type 'sage.rings.integer.Integer'>
            """
            highest_weight = self(highest_weight)
            if not highest_weight.is_dominant():
                raise ValueError("the highest weight must be dominant")
            rho = self.rho()
            pr = self.coroot_lattice().positive_roots()
            from sage.rings.integer import Integer
            n = prod(((rho+highest_weight).scalar(x) for x in pr), Integer(1))
            d = prod((rho.scalar(x) for x in pr), Integer(1))
            return Integer(n/d)

    def weil_representation(self) :
        r"""
        OUTPUT:
        
        - A pair of matrices corresponding to T and S.
        """
        disc_bilinear = lambda a, b: (self._dual_basis * vector(QQ, a.lift())) * self._L * (self._dual_basis * vector(QQ, b.lift()))
        disc_quadratic = lambda a: disc_bilinear(a, a) / ZZ(2)
        
        zeta_order = ZZ(lcm([8, 12, prod(self.invariants())] + map(lambda ed: 2 * ed, self.invariants())))
        K = CyclotomicField(zeta_order); zeta = K.gen()

        R = PolynomialRing(K, 'x'); x = R.gen()
#        sqrt2s = (x**2 - 2).factor()
#        if sqrt2s[0][0][0].complex_embedding().real() > 0 :        
#            sqrt2  = sqrt2s[0][0][0]
#        else : 
#            sqrt2  = sqrt2s[0][1]
        Ldet_rts = (x**2 - prod(self.invariants())).factor()
        if Ldet_rts[0][0][0].complex_embedding().real() > 0 :
            Ldet_rt  = Ldet_rts[0][0][0] 
        else :
            Ldet_rt  = Ldet_rts[0][0][0]
                
        Tmat  = diagonal_matrix( K, [zeta**(zeta_order*disc_quadratic(a)) for a in self] )
        Smat = zeta**(zeta_order / 8 * self._L.nrows()) / Ldet_rt  \
               * matrix( K,  [ [ zeta**ZZ(-zeta_order * disc_bilinear(gamma,delta))
                                 for delta in self ]
                               for gamma in self ])
        
        return (Tmat, Smat)

                def rational_catalan_number(self, p, polynomial=False):
                    r"""
                    Return the ``p``-th rational Catalan number
                    associated to ``self``.

                    It is defined by

                    .. MATH::

                        \prod_{i = 1}^n \frac{p + (p(d_i-1)) \mod h)}{d_i},

                    where `d_1, \ldots, d_n` are the degrees and
                    `h` is the Coxeter number. See [STW2016]_
                    for this formula.

                    INPUT:

                    - ``polynomial`` -- optional boolean (default ``False``)
                      if ``True``, return instead the `q`-analogue as a
                      polynomial in `q`

                    REFERENCES:

                    .. [STW2016] C. Stump, H. Thomas, N. Williams.
                       *Cataland II*, in preparation, 2016.

                    EXAMPLES::

                        sage: W = ColoredPermutations(1,3)
                        sage: [W.rational_catalan_number(p) for p in [5,7,8]]
                        [7, 12, 15]

                        sage: W = ColoredPermutations(2,2)
                        sage: [W.rational_catalan_number(p) for p in [7,9,11]]
                        [10, 15, 21]

                    TESTS::

                        sage: W = ColoredPermutations(1,4)
                        sage: W.rational_catalan_number(3, polynomial=True)
                        q^6 + q^4 + q^3 + q^2 + 1
                    """
                    from sage.arith.all import gcd
                    from sage.combinat.q_analogues import q_int

                    h = self.coxeter_number()
                    if not gcd(h,p) == 1:
                        raise ValueError("parameter p = %s is not coprime to the Coxeter number %s" % (p, h))

                    if polynomial:
                        f = q_int
                    else:
                        f = lambda n: n

                    num = prod(f(p + (p * (deg - 1)) % h)
                               for deg in self.degrees())
                    den = prod(f(deg) for deg in self.degrees())
                    return num // den

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 = Zmod(n).cyclotomic_cosets(q); 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
        sage: codes.DuadicCodeOddPair(GF(q),S1,S2)
        ([11, 6] Cyclic Code over GF(3),
         [11, 6] Cyclic Code over GF(3))

    This is consistent with Theorem 6.1.3 in [HP2003]_.
    """
    from .cyclic_code import CyclicCode
    n = len(S1) + len(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).coefficients(sparse=False)]
    coeffs2 = [_lift2smallest_field(c)[0] for c in (g2+j).coefficients(sparse=False)]
    gg1 = P2(coeffs1)
    gg2 = P2(coeffs2)
    gg1 = gcd(gg1, x**n - 1)
    gg2 = gcd(gg2, x**n - 1)
    C1 = CyclicCode(length = n, generator_pol = gg1)
    C2 = CyclicCode(length = n, generator_pol = gg2)
    return C1,C2

def c2(part, t):
    r"""
    Returns the t-Jack scalar product between ``J(part)`` and ``Q(part)``.

    INPUT:

    - ``self`` -- a Jack basis of the symmetric functions
    - ``part`` -- a partition
    - ``t`` -- an optional parameter (default: uses the parameter `t` from the
      Jack basis)

    OUTPUT:

    - a polynomial in the parameter ``t`` which is equal to the scalar
      product of ``J(part)`` and ``Q(part)``

    EXAMPLES::

        sage: from sage.combinat.sf.jack import c2
        sage: t = QQ['t'].gen()
        sage: [c2(p,t) for p in Partitions(3)]
        [6*t^3, 2*t^3 + t^2, t^3 + 3*t^2 + 2*t]
    """
    return prod(
        [t + t * part.arm_lengths(flat=True)[i] + part.leg_lengths(flat=True)[i] for i in range(sum(part))],
        t.parent().one(),
    )

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.all import prod
    return prod([(np+1)*np**(e-1) \
                     for np,e in [(p.absolute_norm(),e) \
                                  for p,e in N.factor()]])

def c1(part, t):
    r"""
    Returns the `t`-Jack scalar product between ``J(part)`` and ``P(part)``.

    INPUT:

    - ``part`` -- a partition
    - ``t`` -- an optional parameter (default: uses the parameter `t` from the
      Jack basis)

    OUTPUT:

    - a polynomial in the parameter ``t`` which is equal to the scalar
      product of ``J(part)`` and ``P(part)``

    EXAMPLES::

        sage: from sage.combinat.sf.jack import c1
        sage: t = QQ['t'].gen()
        sage: [c1(p,t) for p in Partitions(3)]
        [2*t^2 + 3*t + 1, t + 2, 6]
    """
    return prod(
        [1 + t * part.arm_lengths(flat=True)[i] + part.leg_lengths(flat=True)[i] for i in range(sum(part))],
        t.parent().one(),
    )

                def number_of_reflections_of_full_support(self):
                    r"""
                    Return the number of reflections with full
                    support.

                    EXAMPLES::

                        sage: W = Permutations(4)
                        sage: W.number_of_reflections_of_full_support()
                        1

                        sage: W = ColoredPermutations(1,4)
                        sage: W.number_of_reflections_of_full_support()
                        1

                        sage: W = CoxeterGroup("B3")
                        sage: W.number_of_reflections_of_full_support()
                        3

                        sage: W = ColoredPermutations(3,3)
                        sage: W.number_of_reflections_of_full_support()
                        3
                    """
                    n = self.rank()
                    h = self.coxeter_number()
                    l = self.cardinality()
                    return (n * h * prod(d for d in self.codegrees() if d != 0)) // l

    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 = builtins.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 __getitem__(self, support):
        r"""
        Return the cyclically decreasing element associated with ``support``.

        INPUT:

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

        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)
        if not set(support).issubset(set(index_set)) or support == index_set:
            raise ValueError("the support must be a proper subset of the index set")
        if not support:
            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 _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_ring=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 = []
        ker = self.rays().matrix().integer_kernel().matrix()
        for phases in itertools.product(units, repeat=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 _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 _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 _kohnen_phi(self, a, t) :
     ## We use a modified power of a, namely for each prime factor p of a
     ## we use p**floor(v_p(a)/2). This is compensated for in the routine
     ## of \rho
     a_modif = 1
     for (p,e) in a.factor() :
         a_modif = a_modif * p**(e // 2)
     
     res = 0
     for dsq in filter(lambda d: d.is_square(), a.divisors()) :
         d = isqrt(dsq)
         
         for g_diag in itertools.ifilter( lambda diag: prod(diag) == d,
                               itertools.product(*[d.divisors() for _ in xrange(t.nrows())]) ) :
             for subents in itertools.product(*[xrange(r) for (j,r) in enumerate(g_diag) for _ in xrange(j)]) :
                 columns = [subents[(j * (j - 1)) // 2:(j * (j + 1)) // 2] for j in xrange(t.nrows())]
                 g = diagonal_matrix(list(g_diag))
                 for j in xrange(t.nrows()) :
                     for i in xrange(j) :
                         g[i,j] = columns[j][i]
                      
                 ginv = g.inverse()   
                 tg = ginv.transpose() * t * ginv
                 try :
                     tg= matrix(ZZ, tg)
                 except :
                     continue
                 if any(tg[i,i] % 2 == 1 for i in xrange(tg.nrows())) :
                     continue
                 
                 tg.set_immutable()
                     
                 res = res + self._kohnen_rho(tg, a // dsq)
     
     return a_modif * res

    def selmer_group_iterator(self, S, m, proof=True):
        r"""
        Return an iterator through elements of the finite group `\QQ(S,m)`.

        INPUT:

        - ``S`` -- a set of primes

        - ``m`` -- a positive integer

        - ``proof`` -- ignored

        OUTPUT:

        An iterator yielding the distinct elements of `\QQ(S,m)`.  See
        the docstring for :meth:`selmer_group` for more information.

        EXAMPLES::

            sage: list(QQ.selmer_group_iterator((), 2))
            [1, -1]
            sage: list(QQ.selmer_group_iterator((2,), 2))
            [1, 2, -1, -2]
            sage: list(QQ.selmer_group_iterator((2,3), 2))
            [1, 3, 2, 6, -1, -3, -2, -6]
            sage: list(QQ.selmer_group_iterator((5,), 2))
            [1, 5, -1, -5]
        """
        KSgens, ords = self.selmer_group(S=S, m=m, proof=proof, orders=True)
        one = self.one()
        from sage.misc.all import prod
        from itertools import product
        for ev in product(*[range(o) for o in ords]):
            yield prod((p**e for p,e in zip(KSgens, ev)), one)

def find_unique_orbit_lucas(k, G):
    '''
    Lucas functions (Pell's conic) variant of Rains' algorithm.
    Only works with k = F_{p^n}, p ≠ 2, and G ⊂ ℤ/m*, where

    1. <p> ⊂ ℤ/m* is of order r⋅2,
    2. gcd(n, 2) = 1,
    3. ℤ/m* = <p^2> × G.

    Return a Gaussian period of m-th roots of unity of k.
    '''
    m = G[0][0].parent().order()
    q = k.cardinality()
    assert((q+1) % m == 0)
    cofactor = (q + 1) // m
    sqtest =  (q - 1) // 2
    fact = m.factor()

    # find an m-th root of unity
    zeta = 2
    while any(lucas_pow(zeta, m // f[0]) == 2 for f in fact):
        a = k.random_element()
        if (a**2 - 4)**sqtest == 1:
            continue
        zeta = lucas_pow(a, cofactor)
        
    # return the Gaussian period
    # ... lovely combinatorial iterators!
    return sum(lucas_pow(zeta, prod(g**e for (g, _), e in zip(G, exps)).lift())
               for exps in CProd(map(lambda (_,x): range(x), G)))/2

def q_factorial(n, q=None):
    """
    Returns the `q`-analogue of the factorial `n!`.

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

    EXAMPLES::

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

    The `q`-analogue of `n!` is only defined for `n` a non-negative
    integer (:trac:`11411`)::

        sage: q_factorial(-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(i, q) for i in range(1, n+1)])
    else:
        raise ValueError("Argument (%s) must be a nonnegative integer." %n)

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 _sage_(self):
        """
        EXAMPLES::

            sage: m = lie('[[1,0,3,3],[12,4,-4,7],[-1,9,8,0],[3,-5,-2,9]]') # optional - lie
            sage: m.sage()  # optional - lie
            [ 1  0  3  3]
            [12  4 -4  7]
            [-1  9  8  0]
            [ 3 -5 -2  9]

        """
        t = self.type()
        if t == "grp":
            raise ValueError("cannot convert Lie groups to native Sage objects")
        elif t == "mat":
            import sage.matrix.constructor

            return sage.matrix.constructor.matrix(eval(str(self).replace("\n", "").strip()))
        elif t == "pol":
            from sage.rings.all import PolynomialRing, QQ

            # Figure out the number of variables
            s = str(self)
            open_bracket = s.find("[")
            close_bracket = s.find("]")
            nvars = len(s[open_bracket:close_bracket].split(","))

            # create the polynomial ring
            R = PolynomialRing(QQ, nvars, "x")
            x = R.gens()
            pol = R(0)

            # Split up the polynomials into terms
            terms = []
            for termgrp in s.split(" - "):
                # The first entry in termgrp has
                # a negative coefficient
                termgrp = "-" + termgrp.strip()
                terms += termgrp.split("+")
            # Make sure we don't accidentally add a negative
            # sign to the first monomial
            if s[0] != "-":
                terms[0] = terms[0][1:]

            # go through all the terms in s
            for term in terms:
                xpos = term.find("X")
                coef = eval(term[:xpos].strip())
                exps = eval(term[xpos + 1 :].strip())
                monomial = prod([x[i] ** exps[i] for i in range(nvars)])
                pol += coef * monomial

            return pol
        elif t == "tex":
            return repr(self)
        elif t == "vid":
            return None
        else:
            return ExpectElement._sage_(self)