def _naive_ideal(self, ring):
        r"""
        Return the "naive" subideal.

        INPUT:

        - ``ring`` -- the ambient ring of the ideal.

        OUTPUT:

        A subideal of the toric ideal in the polynomial ring ``ring``.

        EXAMPLES::

            sage: A = matrix([[1,1,1],[0,1,2]])
            sage: IA = ToricIdeal(A)
            sage: IA.ker()
            Free module of degree 3 and rank 1 over Integer Ring
            User basis matrix:
            [ 1 -2  1]
            sage: IA._naive_ideal(IA.ring())
            Ideal (-z1^2 + z0*z2) of Multivariate Polynomial Ring in z0, z1, z2 over Rational Field
        """
        x = ring.gens()
        binomials = []
        for row in self.ker().matrix().rows():
            xpos = prod(x[i]**max( row[i],0) for i in range(0,len(x)))
            xneg = prod(x[i]**max(-row[i],0) for i in range(0,len(x)))
            binomials.append(xpos - xneg)
        return ring.ideal(binomials)

    def all_vv_verifs(self, u, v, dim):
        a,b = set(u).union(set(v))
        a1,b1 = symbolic_max_plus_matrices_band(dim, 2, 'v', 'v', typ='sym')
        d1 = {a:a1, b:b1}
        a2,b2 = symbolic_max_plus_matrices_band(dim, 2, 'v', 'v', typ='full')
        d2 = {a:a2, b:b2}

        mu1 = prod(d1[i] for i in u)
        mu2 = prod(d2[i] for i in u)
        mv1 = prod(d1[i] for i in v)
        mv2 = prod(d2[i] for i in v)

        is_identity = is_vv_identity(u,v,dim)

        self.assertEqual(is_identity,
                         mu1 == mv1,
                         'u = {} v = {} dim = {}'.format(u,v,dim))
        self.assertEqual(mu1 == mv1, mu2 == mv2,
                         'u = {} v = {} dim = {}'.format(u,v,dim))
        if is_identity:
            elts = [random_integer_max_plus_matrices_band(
                       dim, -1000, 1000, ord('v'), ord('v')) \
                    for _ in range(100)]
            zo = {a:0, b:1}
            t0 = tuple(zo[i] for i in u)
            t1 = tuple(zo[i] for i in v)
            self.assertTrue(is_relation(t0,t1, elts, True),
                            'u = {} v = {} dim = {}'.format(u,v,dim))
            self.assertTrue(is_relation(t0,t1,elts,False),
                            'u = {} v = {} dim = {}'.format(u,v,dim))

 def d_vector(self):
     n = self.parent().rk
     one = self.parent().ambient_field()(1)
     factors = self.lift_to_field().factor()
     initial = []
     non_initial = []
     [(initial if x[1] > 0 and len(x[0].monomials()) == 1 else non_initial).append(x[0]**x[1]) for x in factors]
     initial = prod(initial+[one]).numerator()
     non_initial = prod(non_initial+[one]).denominator()
     v1 = vector(non_initial.exponents()[0][:n])
     v2 = vector(initial.exponents()[0][:n])
     return tuple(v1-v2)

        def scalar(self, y):
            r"""
            Return the scalar product of ``self`` and ``y``.

            The scalar product is given by

            .. MATH::

                (I_{\lambda}, I_{\mu}) = \delta_{\lambda,\mu}
                \frac{1}{a_{\lambda}},

            where `a_{\lambda}` is given by

            .. MATH::

                a_{\lambda} = q^{|\lambda| + 2 n(\lambda)} \prod_k
                \prod_{i=1}^{l_k} (1 - q^{-i})

            where `n(\lambda) = \sum_i (i - 1) \lambda_i` and
            `\lambda = (1^{l_1}, 2^{l_2}, \ldots, m^{l_m})`.

            Note that `a_{\lambda}` can be interpreted as the number
            of automorphisms of a certain object in a category
            corresponding to `\lambda`. See Lemma 2.8 in [Schiffmann]_
            for details.

            EXAMPLES::

                sage: R.<q> = ZZ[]
                sage: H = HallAlgebra(R, q)
                sage: H[1].scalar(H[1])
                1/(q - 1)
                sage: H[2].scalar(H[2])
                1/(q^2 - q)
                sage: H[2,1].scalar(H[2,1])
                1/(q^5 - 2*q^4 + q^3)
                sage: H[1,1,1,1].scalar(H[1,1,1,1])
                1/(q^16 - q^15 - q^14 + 2*q^11 - q^8 - q^7 + q^6)
                sage: H.an_element().scalar(H.an_element())
                (4*q^2 + 9)/(q^2 - q)
            """
            q = self.parent()._q
            f = lambda la: ~( q**(sum(la) + 2*la.weighted_size())
                              * prod(prod((1 - q**-i) for i in range(1,k+1))
                                     for k in la.to_exp()) )
            y = self.parent()(y)
            ret = q.parent().zero()
            for mx, cx in self:
                cy = y.coefficient(mx)
                if cy != 0:
                    ret += cx * cy * f(mx)
            return ret

    def defining_polynomials(self):
        """
        Return the pair of products of cyclotomic polynomials.

        EXAMPLES::

            sage: from sage.modular.hypergeometric_motive import HypergeometricData as Hyp
            sage: Hyp(alpha_beta=([1/4,3/4],[0,0])).defining_polynomials()
            (x^2 + 1, x^2 - 2*x + 1)
        """
        up = prod(cyclotomic_polynomial(d) for d in self._cyclo_up)
        down = prod(cyclotomic_polynomial(d) for d in self._cyclo_down)
        return (up, down)

    def quantum_determinant(self, u=None):
        """
        Return the quantum determinant of ``self``.

        The quantum determinant is defined by:

        .. MATH::

            \operatorname{qdet}(u) = \sum_{\sigma \in S_n} (-1)^{\sigma}
            \prod_{k=1}^n T_{\sigma(k),k}(u - k + 1).

        EXAMPLES::

            sage: Y = Yangian(QQ, 2, 2)
            sage: Y.quantum_determinant()
            u^4 + (-2 + t(1)[1,1] + t(1)[2,2])*u^3
             + (1 - t(1)[1,1] + t(1)[1,1]*t(1)[2,2] - t(1)[1,2]*t(1)[2,1]
                - 2*t(1)[2,2] + t(2)[1,1] + t(2)[2,2])*u^2
             + (-t(1)[1,1]*t(1)[2,2] + t(1)[1,1]*t(2)[2,2]
                + t(1)[1,2]*t(1)[2,1] - t(1)[1,2]*t(2)[2,1]
                - t(1)[2,1]*t(2)[1,2] + t(1)[2,2] + t(1)[2,2]*t(2)[1,1]
                - t(2)[1,1] - t(2)[2,2])*u
             - t(1)[1,1]*t(2)[2,2] + t(1)[1,2]*t(2)[2,1] + t(2)[1,1]*t(2)[2,2]
                - t(2)[1,2]*t(2)[2,1] + t(2)[2,2]
        """
        if u is None:
            u = PolynomialRing(self.base_ring(), 'u').gen(0)
        from sage.combinat.permutation import Permutations
        n = self._n
        return sum(p.sign() * prod(self.defining_polynomial(p[k], k+1, u - k)
                                   for k in range(n))
                   for p in Permutations(n))

    def braid(self, n, K=QQ, names=None):
        r"""
        The braid arrangement.

        INPUT:

        - ``n`` -- integer

        - ``K`` -- field (default: ``QQ``)

        - ``names`` -- tuple of strings or ``None`` (default); the
          variable names for the ambient space

        OUTPUT:

        The hyperplane arrangement consisting of the `n(n-1)/2`
        hyperplanes `\{ x_i - x_j = 0 : 1 \leq i \leq j \leq n \}`.

        EXAMPLES::

            sage: hyperplane_arrangements.braid(4)
            Arrangement of 6 hyperplanes of dimension 4 and rank 3
        """
        x = polygen(QQ, 'x')
        A = self.graphical(graphs.CompleteGraph(n), K, names=names)
        charpoly = prod(x-i for i in range(n))
        A.characteristic_polynomial.set_cache(charpoly)
        return A

    def e_one_star_patch(self, v, n):
        r"""
        Return the n-th iterated patch of normal vector v.

        INPUT:

        - ``v`` -- vector, the normal vector
        - ``n`` -- integer

        EXAMPLES::

            sage: from slabbe.mult_cont_frac import FullySubtractive
            sage: FullySubtractive().e_one_star_patch((1,e,pi), 4)
            Patch of 21 faces
        """
        from sage.combinat.e_one_star import E1Star, Patch, Face
        from sage.misc.misc_c import prod
        if v is None:
            v = (random(), random(), random())
        it = self.coding_iterator(v)
        keys = [next(it) for _ in range(n)]
        D = self.dual_substitutions()
        L = prod(D[key] for key in reversed(keys))
        dual_sub = E1Star(L)
        cube = Patch([Face((1,0,0),1), Face((0,1,0),2), Face((0,0,1),3)])
        return dual_sub(cube)

    def tree_factorial(self):
        """
        Returns the tree-factorial of ``self``

        Definition:

        The tree-factorial `T!` of a tree `T` is the product `\prod_{v\in
        T}\#\mbox{children}(v)`.

        EXAMPLES::

            sage: LT = LabelledOrderedTrees()
            sage: t = LT([LT([],label=6),LT([],label=1)],label=9)
            sage: t.tree_factorial()
            3

            sage: BinaryTree([[],[[],[]]]).tree_factorial()
            15

        TESTS::

            sage: BinaryTree().tree_factorial()
            1
        """
        nb = self.node_number()
        if nb <= 1:
            return 1
        else:
            return nb*prod(s.tree_factorial() for s in self)

def m_to_s_stat(R, I, K):
    r"""
    Returns the statistic for the expansion of the Monomial basis element indexed by two
    compositions, as in formula (36) of Tevlin's "Noncommutative Analogs of Monomial Symmetric
    Functions, Cauchy Identity, and Hall Scalar Product".

    INPUT:

    - ``R`` -- A ring
    - ``I``, ``K`` -- compositions

    OUTPUT:

    - An integer

    EXAMPLES::

        sage: from sage.combinat.ncsf_qsym.combinatorics import m_to_s_stat
        sage: m_to_s_stat(QQ,Composition([2,1]), Composition([1,1,1]))
        -1
        sage: m_to_s_stat(QQ,Composition([3]), Composition([1,2]))
        -2
    """
    stat = 0
    for J in Compositions(I.size()):
        if (I.is_finer(J) and  K.is_finer(J)):
            pvec = [0] + Composition(I).refinement_splitting_lengths(J).partial_sums()
            pp = prod( R( len(I) - pvec[i] ) for i in range( len(pvec)-1 ) )
            stat += R((-1)**(len(I)-len(K)) / pp * coeff_lp(K,J))
    return stat

    def bch_to_grs(self):
        r"""
        Returns the underlying GRS code from which ``self`` was derived.

        EXAMPLES::

            sage: C = codes.BCHCode(GF(2), 15, 3)
            sage: RS = C.bch_to_grs()
            sage: RS
            [15, 13, 3] Reed-Solomon Code over GF(16)
            sage: C.generator_matrix() * RS.parity_check_matrix().transpose() == 0
            True
        """
        l = self.jump_size()
        b = self.offset()
        n = self.length()
        designed_distance = self.designed_distance()
        grs_dim = n - designed_distance + 1

        alpha = self.primitive_root()
        alpha_l = alpha ** l
        alpha_b = alpha ** b
        evals = [alpha_l ** i for i in range(n)]
        pcm = [alpha_b ** i for i in range(n)]

        multipliers_product = [1/prod([evals[i] - evals[h] for h in range(n) if h != i]) for i in range(n)]
        column_multipliers = [multipliers_product[i]/pcm[i] for i in range(n)]

        return GeneralizedReedSolomonCode(evals, grs_dim, column_multipliers)

def Polya(G, poids):
    """
    Implantation de la formule d'énumération de Pòlya

    INPUT:

    - ``G`` -- un groupe de permutations d'un ensemble `E`
    - ``poids`` -- la liste des poids `w(c)` des éléments `c` d'un ensemble `F`

    Cette fonction renvoie la série génératrice par poids des
    fonctions de `E` dans `F`:

    .. math:: \sum_{f\in E^F} \prod_{e\in E} w(f(e))

    EXAMPLES:

    On calcule le nombre de colliers bicolores à rotation près::

        sage: Polya(CyclicPermutationGroup(5), [1,1])
        8

    Même chose, raffinée par nombre de perles d'une couleur donnée::

        sage: q = QQ['q'].gen()
        sage: Polya(CyclicPermutationGroup(5), [1,q])
        q^5 + q^4 + 2*q^3 + 2*q^2 + q + 1

    .. TODO:: Rajouter ici les autres exemples!

    """
    return sum(prod( p(k, poids) for k in type_cyclique(sigma))
               for sigma in G) / G.cardinality()

        def sum_of_partitions(self, la):
            r"""
            Return the sum over all sets partitions whose shape is ``la``,
            scaled by `\prod_i m_i!` where `m_i` is the multiplicity
            of `i` in ``la``.

            INPUT:

            - ``la`` -- an integer partition

            OUTPUT:

            - an element of ``self``

            EXAMPLES::

                sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                sage: w.sum_of_partitions([2,1,1])
                2*w{{1}, {2}, {3, 4}} + 2*w{{1}, {2, 3}, {4}} + 2*w{{1}, {2, 4}, {3}}
                 + 2*w{{1, 2}, {3}, {4}} + 2*w{{1, 3}, {2}, {4}} + 2*w{{1, 4}, {2}, {3}}
            """
            la = Partition(la)
            c = prod([factorial(_) for _ in la.to_exp()])
            P = SetPartitions()
            return self.sum_of_terms([(P(m), c) for m in SetPartitions(sum(la), la)], distinct=True)

def pure_tensor(P, vects, indices=None):
    r"""
    Return the product of a collection of linear forms corresponding to vectors

    INPUT:

    - ``P`` -- a polynomial ring
    - ``vects`` -- a list of coefficient vectors
    - ``indices`` -- (default: ``None``) a collection of indices for the list
      ``vects`` corresponding to the linear forms which will be included in the
      product.  Indices may be repeated.  If ``None``, then each vector will be
      included in the product once.

    OUTPUT:

    - the corresponding product of linear forms in the polynomial ring ``P``

    EXAMPLES:

        sage: P.<x, y> = PolynomialRing(QQ)
        sage: vects = [[1,0],[0,1],[1,1]]
        sage: pure_tensor(P, vects).factor()
        y * x * (x + y)
        sage: pure_tensor(P, vects, [0, 2, 2]).factor()
        x * (x + y)^2
    """
    if indices is None:
        indices = range(len(vects))
    terms = [linear_form(P, vects[i]) for i in indices]
    return prod(terms, P.one())

def homogenous_symmetric_function(j, x):
    r"""
    Return a complete homogeneous symmetric polynomial
    (:wikipedia:`Complete_homogeneous_symmetric_polynomial`).

    INPUT:

    - ``j`` -- the degree as a nonnegative integer

    - ``x`` -- an iterable of variables

    OUTPUT:

    A polynomial of the common parent of all entries of ``x``

    EXAMPLES::

        sage: from sage.rings.polynomial.omega import homogenous_symmetric_function
        sage: P = PolynomialRing(ZZ, 'X', 3)
        sage: homogenous_symmetric_function(0, P.gens())
        1
        sage: homogenous_symmetric_function(1, P.gens())
        X0 + X1 + X2
        sage: homogenous_symmetric_function(2, P.gens())
        X0^2 + X0*X1 + X1^2 + X0*X2 + X1*X2 + X2^2
        sage: homogenous_symmetric_function(3, P.gens())
        X0^3 + X0^2*X1 + X0*X1^2 + X1^3 + X0^2*X2 +
        X0*X1*X2 + X1^2*X2 + X0*X2^2 + X1*X2^2 + X2^3
    """
    from sage.combinat.integer_vector import IntegerVectors
    from sage.misc.misc_c import prod

    return sum(prod(xx**pp for xx, pp in zip(x, p))
               for p in IntegerVectors(j, length=len(x)))

    def _convert_factors_(self, factors):
        r"""
        Helper method. Try to convert some ``factors`` to an
        element of one of the Cartesian factors and return the product of
        all these factors.

        INPUT:

        - ``factors`` -- a tuple or other iterable.

        OUTPUT:

        An element of this Cartesian product.

        EXAMPLES::

            sage: from sage.rings.asymptotic.growth_group import GrowthGroup
            sage: G = GrowthGroup('x^ZZ * log(x)^QQ * y^QQ')
            sage: e1 = G._convert_factors_([x^2])
            sage: (e1, e1.parent())
            (x^2, Growth Group x^ZZ * log(x)^QQ * y^QQ)
        """
        from sage.misc.misc_c import prod

        def get_factor(data):
            for factor in self.cartesian_factors():
                try:
                    return factor, factor(data)
                except (ValueError, TypeError):
                    pass
            raise ValueError("%s is not in any of the factors of %s" % (data, self))

        return prod(self.cartesian_injection(*get_factor(f)) for f in factors)

    def _homogeneous_generators_noncommutative_variables_zero_Hecke(self, r):
        r"""
        Returns the ``r^{th}`` homogeneous generator, viewed as an element inside the
        affine zero Hecke algebra. This is the sum of all cyclicly decreasing elements
        of order ``r``.

        INPUT:

        - ``r`` -- A positive integer

        OUTPUT:

        - An element of the affine zero Hecke algebra.

        EXAMPLES::

            sage: g = SymmetricFunctions(QQ).kBoundedSubspace(3,1).K_kschur()
            sage: g._homogeneous_generators_noncommutative_variables_zero_Hecke(2)
            T1*T0 + T2*T0 + T0*T3 + T3*T2 + T3*T1 + T2*T1
            sage: g._homogeneous_generators_noncommutative_variables_zero_Hecke(0)
            1
        """
        from sage.combinat.root_system.weyl_group import WeylGroup
        from sage.algebras.iwahori_hecke_algebra import IwahoriHeckeAlgebraT
        W = WeylGroup(['A',self.k,1])
        H = IwahoriHeckeAlgebraT(W, 0, base_ring = self.base_ring())
        Hgens = H.algebra_generators()
        S = [w.reduced_word() for w in W.pieri_factors() if w.length() == r]
        return sum( (prod((Hgens[i] for i in w), 1) for w in S), 0 )

    def Ish(self, n, K=QQ, names=None):
        r"""
        Return the Ish arrangement.

        INPUT:

        - ``n`` -- integer

        - ``K`` -- field (default:``QQ``)

        - ``names`` -- tuple of strings or ``None`` (default); the
          variable names for the ambient space

        OUTPUT:

        The Ish arrangement, which is the set of `n(n-1)` hyperplanes.

        .. MATH::

            \{ x_i - x_j = 0 : 1 \leq i \leq j \leq n \} 
            \cup 
            \{ x_1 - x_j = i : 1 \leq i \leq j \leq n \}.

        EXAMPLES::

            sage: a = hyperplane_arrangements.Ish(3);  a
            Arrangement of 6 hyperplanes of dimension 3 and rank 2
            sage: a.characteristic_polynomial()
            x^3 - 6*x^2 + 9*x
            sage: b = hyperplane_arrangements.Shi(3)
            sage: b.characteristic_polynomial()
            x^3 - 6*x^2 + 9*x

        TESTS::

            sage: a.characteristic_polynomial.clear_cache()  # long time
            sage: a.characteristic_polynomial()              # long time
            x^3 - 6*x^2 + 9*x

        REFERENCES:

        ..  [AR] D. Armstrong, B. Rhoades
            "The Shi arrangement and the Ish arrangement"
            :arxiv:`1009.1655`
        """
        H = make_parent(K, n, names)
        x = H.gens()
        hyperplanes = []
        for i in range(n):
            for j in range(i+1, n):
                hyperplanes.append(x[i] - x[j])
                hyperplanes.append(x[0] - x[j] - (i+1))
        A = H(*hyperplanes)
        x = polygen(QQ, 'x')
        charpoly = x * sum([(-1)**k * stirling_number2(n, n-k) *
                            prod([(x - 1 - j) for j in range(k, n-1)]) for k in range(0, n)])
        A.characteristic_polynomial.set_cache(charpoly)
        return A

            def is_symmetric(self):
                r"""
                Determine if a `NCSym^*` function, expressed in the
                `\mathbf{w}` basis, is symmetric.

                A function `f` in the `\mathbf{w}` basis is a symmetric
                function if it is in the image of `\chi^*`. That is to say we
                have

                .. MATH::

                    f = \sum_{\lambda} c_{\lambda} \prod_i m_i(\lambda)!
                    \sum_{\lambda(A) = \lambda} \mathbf{w}_A

                where the second sum is over all set partitions `A` whose
                shape `\lambda(A)` is equal to `\lambda` and `m_i(\mu)` is
                the multiplicity of `i` in the partition `\mu`.

                OUTPUT:

                - ``True`` if `\lambda(A)=\lambda(B)` implies the coefficients of
                  `\mathbf{w}_A` and `\mathbf{w}_B` are equal, ``False`` otherwise

                EXAMPLES::

                    sage: w = SymmetricFunctionsNonCommutingVariables(QQ).dual().w()
                    sage: elt = w.sum_of_partitions([2,1,1])
                    sage: elt.is_symmetric()
                    True
                    sage: elt -= 3*w.sum_of_partitions([1,1])
                    sage: elt.is_symmetric()
                    True
                    sage: w = SymmetricFunctionsNonCommutingVariables(ZZ).dual().w()
                    sage: elt = w.sum_of_partitions([2,1,1]) / 2
                    sage: elt.is_symmetric()
                    False
                    sage: elt = w[[1,3],[2]]
                    sage: elt.is_symmetric()
                    False
                    sage: elt = w[[1],[2,3]] + w[[1,2],[3]] + 2*w[[1,3],[2]]
                    sage: elt.is_symmetric()
                    False
                """
                d = {}
                R = self.base_ring()
                for A, coeff in self:
                    la = A.shape()
                    exp = prod([factorial(_) for _ in la.to_exp()])
                    if la not in d:
                        if coeff / exp not in R:
                            return False
                        d[la] = [coeff, 1]
                    else:
                        if d[la][0] != coeff:
                            return False
                        d[la][1] += 1
                # Make sure we've seen each set partition of the shape
                return all(d[la][1] == SetPartitions(la.size(), la).cardinality() for la in d)

    def __init__(self, params, asym=False):
        # set parameters
        (self.alpha, self.beta, self.rho, self.rho_f, self.eta, self.bound, self.n, self.k) = params

        self.x0 = ZZ(1)
        
        self.primes = [random_prime(2**self.eta, lbound = 2**(self.eta - 1), proof=False) for i in range(self.n)]
        primes = self.primes
        
        self.x0 = prod(primes)

        # generate CRT coefficients
        self.coeff = [ZZ((self.x0/p_i) * ZZ(Zmod(p_i)(self.x0/p_i)**(-1))) for p_i in primes]

        # generate secret g_i
        self.g = [random_prime(2**self.alpha, proof=False) for i in range(self.n)]

        # generate zs and zs^(-1)
        z = []
        zinv = []

        # generate z and z^(-1)
        if not asym:
            while True:
                z = ZZ.random_element(self.x0)  
                try:
                    zinv = ZZ(Zmod(self.x0)(z)**(-1))
                    break
                except ZeroDivisionError:
                    ''' Error occurred, retry sampling '''

            z, self.zinv = zip(*[(z,zinv) for i in range(self.k)])

        else: # asymmetric version
            for i in range(self.k):
                while True:
                    z_i = ZZ.random_element(self.x0)  
                    try:
                        zinv_i = ZZ(Zmod(self.x0)(z_i)**(-1))
                        break
                    except ZeroDivisionError:
                        ''' Error occurred, retry sampling '''

                z.append(z_i)
                zinv.append(zinv_i)

            self.zinv = zinv

        # generate p_zt
        zk = Zmod(self.x0)(1)
        self.p_zt = 0
        for z_i in z:
            zk *= Zmod(self.x0)(z_i)
        for i in range(self.n):
            self.p_zt += Zmod(self.x0)(ZZ(Zmod(self.primes[i])(self.g[i])**(-1) * Zmod(self.primes[i])(zk)) * ZZ.random_element(2**self.beta) * (self.x0/self.primes[i]))

        self.p_zt = Zmod(self.x0)(self.p_zt)