def elias_upper_bound(n,q,d,algorithm=None):
    r"""
    Returns the Elias upper bound.

    Returns the Elias upper bound for number of elements in the largest
    code of minimum distance `d` in `\GF{q}^n`, cf. [HP2003]_.
    If the method is "gap", it wraps GAP's ``UpperBoundElias``. 

    EXAMPLES::

        sage: codes.bounds.elias_upper_bound(10,2,3)
        232
        sage: codes.bounds.elias_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package)
        232
    """
    _check_n_q_d(n, q, d, field_based=False)
    r = 1-1/q
    if algorithm=="gap":
        gap.load_package("guava")
        ans=gap.eval("UpperBoundElias(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        def ff(n,d,w,q):
            return r*n*d*q**n/((w**2-2*r*n*w+r*n*d)*volume_hamming(n,q,w));
    def get_list(n,d,q):
        I = []
        for i in range(1,int(r*n)+1):
            if i**2-2*r*n*i+r*n*d>0:
                I.append(i)
        return I
    I = get_list(n,d,q)
    bnd = min([ff(n,d,w,q) for w in I])
    return int(bnd)

def plotkin_upper_bound(n,q,d, algorithm=None):
    r"""
    Returns Plotkin upper bound for number of elements in the largest
    code of minimum distance d in `\GF{q}^n`. 

    The algorithm="gap" option wraps Guava's UpperBoundPlotkin.
    
    EXAMPLES::
        sage: plotkin_upper_bound(10,2,3) 
        192
        sage: plotkin_upper_bound(10,2,3,algorithm="gap")  # requires optional GAP package Guava
        192
    """
    if algorithm=="gap":
        ans=gap.eval("UpperBoundPlotkin(%s,%s,%s)"%(n,d,q))
        #print "calling Guava ..."
        return QQ(ans)
    else:
        t = 1 - 1/q
        if (q==2) and (n == 2*d) and (d%2 == 0):
            return 4*d
        elif (q==2) and (n == 2*d + 1) and (d%2 == 1):
            return 4*d + 4
        elif d > t*n:
            return int(d/( d - t*n))
        elif d < t*n + 1:
            fact = (d-1) / t
            if RR(fact)==RR(int(fact)):
                fact = int(fact) + 1
            return int(d/( d - t * fact)) * q**(n - fact)

def elias_upper_bound(n,q,d,algorithm=None):
    r"""
    Returns the Elias upper bound for number of elements in the largest
    code of minimum distance d in `\GF{q}^n`. Wraps
    GAP's UpperBoundElias.
    
    EXAMPLES::
    
        sage: elias_upper_bound(10,2,3)
        232
        sage: elias_upper_bound(10,2,3,algorithm="gap")  # requires optional GAP package Guava
        232

    """
    r = 1-1/q
    if algorithm=="gap":
        #print "calling Guava ..."
        ans=gap.eval("UpperBoundElias(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        def ff(n,d,w,q):
            return r*n*d*q**n/((w**2-2*r*n*w+r*n*d)*volume_hamming(n,q,w));
    def get_list(n,d,q):
        I = []
        for i in range(1,int(r*n)+1):
            if i**2-2*r*n*i+r*n*d>0:
                I.append(i)
            return I
    I = get_list(n,d,q)
    bnd = min([ff(n,d,w,q) for w in I])
    return int(bnd)

    def cpi(self, p):
        """
        Returns the centrally primitive idempotent for the symmetric group
        of order n for the irreducible corresponding indexed by the
        partition p.
        
        EXAMPLES::
        
            sage: QS3 = SymmetricGroupAlgebra(QQ,3)
            sage: QS3.cpi([2,1])
            2/3*[1, 2, 3] - 1/3*[2, 3, 1] - 1/3*[3, 1, 2]
            sage: QS3.cpi([3])
            1/6*[1, 2, 3] + 1/6*[1, 3, 2] + 1/6*[2, 1, 3] + 1/6*[2, 3, 1] + 1/6*[3, 1, 2] + 1/6*[3, 2, 1]
            sage: QS3.cpi([1,1,1])
            1/6*[1, 2, 3] - 1/6*[1, 3, 2] - 1/6*[2, 1, 3] + 1/6*[2, 3, 1] + 1/6*[3, 1, 2] - 1/6*[3, 2, 1]
        """
        if p not in partition.Partitions_n(self.n):
            raise TypeError, "p must be a partition of %s"%self.n

        character_table = eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%self.n))

        cpi = self(0)

        np = partition.Partitions_n(self.n).list()
        np.reverse()
        p_index = np.index(p)
        
        big_coeff = character_table[p_index][0]/factorial(self.n)

        for g in permutation.StandardPermutations_n(self.n):
            cpi += big_coeff * character_table[p_index][np.index(g.inverse().cycle_type())] * self(g)

        return cpi

def codesize_upper_bound(n,d,q,algorithm=None):
    r"""
    This computes the minimum value of the upper bound using the 
    algorithms of Singleton, Hamming, Plotkin and Elias.

    If algorithm="gap" then this returns the best known upper 
    bound `A(n,d)=A_q(n,d)` for the size of a code of length n, 
    minimum distance d over a field of size q. The function first 
    checks for trivial cases (like d=1 or n=d), and if the value 
    is in the built-in table. Then it calculates the minimum value 
    of the upper bound using the algorithms of Singleton, Hamming, 
    Johnson, Plotkin and Elias. If the code is binary, 
    `A(n, 2\ell-1) = A(n+1,2\ell)`, so the function
    takes the minimum of the values obtained from all algorithms for the
    parameters `(n, 2\ell-1)` and `(n+1, 2\ell)`. This
    wraps GUAVA's UpperBound( n, d, q ).
    
    EXAMPLES::
        sage: codesize_upper_bound(10,3,2)
        93
        sage: codesize_upper_bound(10,3,2,algorithm="gap")  # requires optional GAP package Guava
        85

    """
    if algorithm=="gap":
        return int(gap.eval("UpperBound(%s,%s,%s)"%( n, d, q )))
    else:
        eub = elias_upper_bound(n,q,d)
        gub = griesmer_upper_bound(n,q,d)
        hub = hamming_upper_bound(n,q,d)
        pub = plotkin_upper_bound(n,q,d)
        sub = singleton_upper_bound(n,q,d)
        return min([eub,gub,hub,pub,sub])

def griesmer_upper_bound(n,q,d,algorithm=None):
    r"""
    Returns the Griesmer upper bound for number of elements in the
    largest code of minimum distance d in `\GF{q}^n`.
    Wraps GAP's UpperBoundGriesmer.
    
    EXAMPLES::
        sage: griesmer_upper_bound(10,2,3)
        128
        sage: griesmer_upper_bound(10,2,3,algorithm="gap")  # requires optional GAP package Guava
        128
    """
    if algorithm=="gap":
        #print "calling Guava ..."
        ans=gap.eval("UpperBoundGriesmer(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        den = 1
        s = 0
        k = 0
        add = 0
        while s <= n:
            if not(add == 1):
                if d%den==0:
                    add = int(d/den)
                else: 
                    add = int(d/den)+1
            s = s + add
            den = den * q
            k = k + 1
        return q**(k-1)

def RandomLinearCodeGuava(n,k,F):
    r"""
    The method used is to first construct a `k \times n` matrix of the block
    form `(I,A)`, where `I` is a `k \times k` identity matrix and `A` is a
    `k \times (n-k)` matrix constructed using random elements of `F`. Then
    the columns are permuted using a randomly selected element of the symmetric
    group `S_n`.

    INPUT:
        Integers `n,k`, with `n>k>1`.

    OUTPUT:
        Returns a "random" linear code with length n, dimension k over field F.

    EXAMPLES::
        sage: C = RandomLinearCodeGuava(30,15,GF(2)); C      # optional - gap_packages (Guava package)
        Linear code of length 30, dimension 15 over Finite Field of size 2
        sage: C = RandomLinearCodeGuava(10,5,GF(4,'a')); C      # optional - gap_packages (Guava package)
        Linear code of length 10, dimension 5 over Finite Field in a of size 2^2

    AUTHOR: David Joyner (11-2005)
    """
    current_randstate().set_seed_gap()

    q = F.order()
    gap.eval("C:=RandomLinearCode("+str(n)+","+str(k)+", GF("+str(q)+"))")
    gap.eval("G:=GeneratorMat(C)")
    k = int(gap.eval("Length(G)"))
    n = int(gap.eval("Length(G[1])"))
    G = [[gfq_gap_to_sage(gap.eval("G[%s][%s]" % (i,j)),F) for j in range(1,n+1)] for i in range(1,k+1)]
    MS = MatrixSpace(F,k,n)
    return LinearCode(MS(G))

    def cardinality(self):
        """
        EXAMPLES::

            sage: S = [1,2,3,4,5]
            sage: UnorderedTuples(S,2).cardinality()
            15
        """
        ans=gap.eval("NrUnorderedTuples(%s,%s)"%(self._index_list,ZZ(self.k)))
        return ZZ(ans)

def griesmer_upper_bound(n,q,d,algorithm=None):
    r"""
    Returns the Griesmer upper bound.

    Returns the Griesmer upper bound for the number of elements in a
    largest linear code of minimum distance `d` in `\GF{q}^n`, cf. [HP2003]_.
    If the method is "gap", it wraps GAP's ``UpperBoundGriesmer``. 

    The bound states:

    .. MATH::

        `n\geq \sum_{i=0}^{k-1} \lceil d/q^i \rceil.`


    EXAMPLES:

    The bound is reached for the ternary Golay codes::

        sage: codes.bounds.griesmer_upper_bound(12,3,6)
        729
        sage: codes.bounds.griesmer_upper_bound(11,3,5)
        729

    ::

        sage: codes.bounds.griesmer_upper_bound(10,2,3)
        128
        sage: codes.bounds.griesmer_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package)
        128

    TESTS::

        sage: codes.bounds.griesmer_upper_bound(11,3,6)
        243
        sage: codes.bounds.griesmer_upper_bound(11,3,6)
        243
    """
    _check_n_q_d(n, q, d)
    if algorithm=="gap":
        gap.load_package("guava")
        ans=gap.eval("UpperBoundGriesmer(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        #To compute the bound, we keep summing up the terms on the RHS
        #until we start violating the inequality.
        from sage.functions.other import ceil
        den = 1
        s = 0
        k = 0
        while s <= n:
            s += ceil(d/den)
            den *= q
            k = k + 1
        return q**(k-1)

    def list(self):
        """
        EXAMPLES::

            sage: S = [1,2]
            sage: UnorderedTuples(S,3).list()
            [[1, 1, 1], [1, 1, 2], [1, 2, 2], [2, 2, 2]]
            sage: UnorderedTuples(["a","b","c"],2).list()
            [['a', 'a'], ['a', 'b'], ['a', 'c'], ['b', 'b'], ['b', 'c'], ['c', 'c']]
        """
        ans=gap.eval("UnorderedTuples(%s,%s)"%(self._index_list,ZZ(self.k)))
        return [[self.S[i] for i in l] for l in eval(ans)]

    def list(self):
        """
        EXAMPLES::

            sage: S = [1,2]
            sage: UnorderedTuples(S,3).list()
            [[1, 1, 1], [1, 1, 2], [1, 2, 2], [2, 2, 2]]
            sage: UnorderedTuples(["a","b","c"],2).list()
            [['a', 'a'], ['a', 'b'], ['a', 'c'], ['b', 'b'], ['b', 'c'], ['c', 'c']]
        """
        ans=gap.eval("UnorderedTuples(%s,%s)"%(self._index_list,ZZ(self.k)))
        return map(lambda l: map(lambda i: self.S[i], l), eval(ans))

    def cardinality(self):
        """
        Returns the size of combinations(mset,k). IMPLEMENTATION: Wraps
        GAP's NrCombinations.

        EXAMPLES::

            sage: mset = [1,1,2,3,4,4,5]
            sage: Combinations(mset,2).cardinality()
            12
        """
        items = map(self.mset.index, self.mset)
        return ZZ(gap.eval("NrCombinations({},{})".format( items,ZZ(self.k))) )

def codesize_upper_bound(n,d,q,algorithm=None):
    r"""
    This computes the minimum value of the upper bound using the
    methods of Singleton, Hamming, Plotkin, and Elias.

    If algorithm="gap" then this returns the best known upper
    bound `A(n,d)=A_q(n,d)` for the size of a code of length n,
    minimum distance d over a field of size q. The function first
    checks for trivial cases (like d=1 or n=d), and if the value
    is in the built-in table. Then it calculates the minimum value
    of the upper bound using the algorithms of Singleton, Hamming,
    Johnson, Plotkin and Elias. If the code is binary,
    `A(n, 2\ell-1) = A(n+1,2\ell)`, so the function
    takes the minimum of the values obtained from all algorithms for the
    parameters `(n, 2\ell-1)` and `(n+1, 2\ell)`. This
    wraps GUAVA's (i.e. GAP's package Guava) UpperBound( n, d, q ).

    If algorithm="LP" then this returns the Delsarte (a.k.a. Linear
    Programming) upper bound.

    EXAMPLES::

        sage: codes.bounds.codesize_upper_bound(10,3,2)
        93
        sage: codes.bounds.codesize_upper_bound(24,8,2,algorithm="LP")
        4096
        sage: codes.bounds.codesize_upper_bound(10,3,2,algorithm="gap")  # optional - gap_packages (Guava package)
        85
        sage: codes.bounds.codesize_upper_bound(11,3,4,algorithm=None)
        123361
        sage: codes.bounds.codesize_upper_bound(11,3,4,algorithm="gap")  # optional - gap_packages (Guava package)
        123361
        sage: codes.bounds.codesize_upper_bound(11,3,4,algorithm="LP")
        109226

    """
    if algorithm=="gap":
        gap.load_package('guava')
        return int(gap.eval("UpperBound(%s,%s,%s)"%( n, d, q )))
    if algorithm=="LP":
        return int(delsarte_bound_hamming_space(n,d,q))
    else:
        eub = elias_upper_bound(n,q,d)
        gub = griesmer_upper_bound(n,q,d)
        hub = hamming_upper_bound(n,q,d)
        pub = plotkin_upper_bound(n,q,d)
        sub = singleton_upper_bound(n,q,d)
        return min([eub,gub,hub,pub,sub])

def BinaryReedMullerCode(r,k):
    r"""
    The binary 'Reed-Muller code' with dimension k and
    order r is a code with length `2^k` and minimum distance `2^k-r`
    (see for example, section 1.10 in [HP]_). By definition, the
    `r^{th}` order binary Reed-Muller code of length `n=2^m`, for
    `0 \leq r \leq m`, is the set of all vectors `(f(p)\ |\ p \\in GF(2)^m)`,
    where `f` is a multivariate polynomial of degree at most `r`
    in `m` variables.

    INPUT:

    - ``r, k`` -- positive integers with `2^k>r`.

    OUTPUT:

    Returns the binary 'Reed-Muller code' with dimension `k` and order `r`.

    EXAMPLE::

        sage: C = codes.BinaryReedMullerCode(2,4); C  # optional - gap_packages (Guava package)
        Linear code of length 16, dimension 11 over Finite Field of size 2
        sage: C.minimum_distance()              # optional - gap_packages (Guava package)
        4
        sage: C.gen_mat()                       # optional - gap_packages (Guava package)
        [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
        [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1]
        [0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1]
        [0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1]
        [0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]
        [0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1]
        [0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 1]
        [0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1]
        [0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1]
        [0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1]
        [0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1]

    AUTHOR: David Joyner (11-2005)
    """
    F = GF(2)
    gap.load_package("guava")
    gap.eval("C:=ReedMullerCode("+str(r)+", "+str(k)+")")
    gap.eval("G:=GeneratorMat(C)")
    k = int(gap.eval("Length(G)"))
    n = int(gap.eval("Length(G[1])"))
    G = [[gfq_gap_to_sage(gap.eval("G["+str(i)+"]["+str(j)+"]"),F) for j in range(1,n+1)] for i in range(1,k+1)]
    MS = MatrixSpace(F,k,n)
    return LinearCode(MS(G))

    def cpi(self, p):
        """
        Return the centrally primitive idempotent for the symmetric group
        of order `n` corresponding to the irreducible representation
        indexed by the partition ``p``.

        EXAMPLES::

            sage: QS3 = SymmetricGroupAlgebra(QQ,3)
            sage: QS3.cpi([2,1])
            2/3*[1, 2, 3] - 1/3*[2, 3, 1] - 1/3*[3, 1, 2]
            sage: QS3.cpi([3])
            1/6*[1, 2, 3] + 1/6*[1, 3, 2] + 1/6*[2, 1, 3] + 1/6*[2, 3, 1] + 1/6*[3, 1, 2] + 1/6*[3, 2, 1]
            sage: QS3.cpi([1,1,1])
            1/6*[1, 2, 3] - 1/6*[1, 3, 2] - 1/6*[2, 1, 3] + 1/6*[2, 3, 1] + 1/6*[3, 1, 2] - 1/6*[3, 2, 1]

            sage: QS0 = SymmetricGroupAlgebra(QQ, 0)
            sage: QS0.cpi(Partition([]))
            []

        TESTS::

            sage: QS3.cpi([2,2])
            Traceback (most recent call last):
            ...
            TypeError: p (= [2, 2]) must be a partition of n (= 3)
        """
        if p not in partition.Partitions_n(self.n):
            raise TypeError("p (= {p}) must be a partition of n (= {n})".format(p=p, n=self.n))

        character_table = eval(gap.eval("Display(Irr(SymmetricGroup(%d)));"%self.n))

        cpi = self.zero()

        np = partition.Partitions_n(self.n).list()
        np.reverse()
        p_index = np.index(p)

        big_coeff = character_table[p_index][0] / factorial(self.n)

        character_row = character_table[p_index]
        dct = { g : big_coeff * character_row[np.index(g.cycle_type())]
                for g in permutation.StandardPermutations_n(self.n) }

        return self._from_dict(dct)

def QuasiQuadraticResidueCode(p):
    r"""
    A (binary) quasi-quadratic residue code (or QQR code), as defined by
    Proposition 2.2 in [BM]_, has a generator matrix in the block form `G=(Q,N)`.
    Here `Q` is a `p \times p` circulant matrix whose top row
    is `(0,x_1,...,x_{p-1})`, where `x_i=1` if and only if `i`
    is a quadratic residue `\mod p`, and `N` is a `p \times p` circulant
    matrix whose top row is `(0,y_1,...,y_{p-1})`, where `x_i+y_i=1` for all `i`.

    INPUT:

    - ``p`` -- a prime `>2`.

    OUTPUT:

    Returns a QQR code of length `2p`.

    EXAMPLES::

        sage: C = codes.QuasiQuadraticResidueCode(11); C   # optional - gap_packages (Guava package)
        Linear code of length 22, dimension 11 over Finite Field of size 2

    REFERENCES:

    .. [BM] Bazzi and Mitter, {\it Some constructions of codes from group actions}, (preprint
      March 2003, available on Mitter's MIT website).

    .. [Jresidue] \D. Joyner, {\it On quadratic residue codes and hyperelliptic curves},
      (preprint 2006)

    These are self-orthogonal in general and self-dual when $p \\equiv 3 \\pmod 4$.

    AUTHOR: David Joyner (11-2005)
    """
    F = GF(2)
    gap.load_package("guava")
    gap.eval("C:=QQRCode(" + str(p) + ")")
    gap.eval("G:=GeneratorMat(C)")
    k = int(gap.eval("Length(G)"))
    n = int(gap.eval("Length(G[1])"))
    G = [[gfq_gap_to_sage(gap.eval("G[%s][%s]" % (i, j)), F)
          for j in range(1, n + 1)] for i in range(1, k + 1)]
    MS = MatrixSpace(F, k, n)
    return LinearCode(MS(G))

def QuasiQuadraticResidueCode(p):
    r"""
    A (binary) quasi-quadratic residue code (or QQR code).

    Follows the definition of Proposition 2.2 in [BM]. The code has a generator
    matrix in the block form `G=(Q,N)`. Here `Q` is a `p \times p` circulant
    matrix whose top row is `(0,x_1,...,x_{p-1})`, where `x_i=1` if and only if
    `i` is a quadratic residue `\mod p`, and `N` is a `p \times p` circulant
    matrix whose top row is `(0,y_1,...,y_{p-1})`, where `x_i+y_i=1` for all
    `i`.

    INPUT:

    - ``p`` -- a prime `>2`.

    OUTPUT:

    Returns a QQR code of length `2p`.

    EXAMPLES::

        sage: C = codes.QuasiQuadraticResidueCode(11); C   # optional - gap_packages (Guava package)
        [22, 11] linear code over GF(2)

    These are self-orthogonal in general and self-dual when $p \\equiv 3 \\pmod 4$.

    AUTHOR: David Joyner (11-2005)
    """
    if not is_package_installed('gap_packages'):
        raise PackageNotFoundError('gap_packages')
    F = GF(2)
    gap.load_package("guava")
    gap.eval("C:=QQRCode(" + str(p) + ")")
    gap.eval("G:=GeneratorMat(C)")
    k = int(gap.eval("Length(G)"))
    n = int(gap.eval("Length(G[1])"))
    G = [[gfq_gap_to_sage(gap.eval("G[%s][%s]" % (i, j)), F)
          for j in range(1, n + 1)] for i in range(1, k + 1)]
    MS = MatrixSpace(F, k, n)
    return LinearCode(MS(G))

def RandomLinearCodeGuava(n, k, F):
    r"""
    The method used is to first construct a `k \times n` matrix of the block
    form `(I,A)`, where `I` is a `k \times k` identity matrix and `A` is a
    `k \times (n-k)` matrix constructed using random elements of `F`. Then
    the columns are permuted using a randomly selected element of the symmetric
    group `S_n`.

    INPUT:

    - ``n,k`` -- integers with `n>k>1`.

    OUTPUT:

    Returns a "random" linear code with length `n`, dimension `k` over field `F`.

    EXAMPLES::

        sage: C = codes.RandomLinearCodeGuava(30,15,GF(2)); C      # optional - gap_packages (Guava package)
        [30, 15] linear code over GF(2)
        sage: C = codes.RandomLinearCodeGuava(10,5,GF(4,'a')); C      # optional - gap_packages (Guava package)
        [10, 5] linear code over GF(4)

    AUTHOR: David Joyner (11-2005)
    """
    current_randstate().set_seed_gap()

    q = F.order()
    if not is_package_installed('gap_packages'):
        raise PackageNotFoundError('gap_packages')
    gap.load_package("guava")
    gap.eval("C:=RandomLinearCode("+str(n)+","+str(k)+", GF("+str(q)+"))")
    gap.eval("G:=GeneratorMat(C)")
    k = int(gap.eval("Length(G)"))
    n = int(gap.eval("Length(G[1])"))
    G = [[gfq_gap_to_sage(gap.eval("G[%s][%s]" % (i, j)), F)
          for j in range(1, n + 1)] for i in range(1, k + 1)]
    MS = MatrixSpace(F, k, n)
    return LinearCode(MS(G))

def plotkin_upper_bound(n,q,d, algorithm=None):
    r"""
    Returns the Plotkin upper bound.

    Returns the Plotkin upper bound for the number of elements in a largest
    code of minimum distance `d` in `\GF{q}^n`.
    More precisely this is a generalization of Plotkin's result for `q=2`
    to bigger `q` due to Berlekamp.

    The ``algorithm="gap"`` option wraps Guava's ``UpperBoundPlotkin``.

    EXAMPLES::

        sage: codes.bounds.plotkin_upper_bound(10,2,3)
        192
        sage: codes.bounds.plotkin_upper_bound(10,2,3,algorithm="gap")  # optional - gap_packages (Guava package)
        192
    """
    _check_n_q_d(n, q, d, field_based=False)
    if algorithm=="gap":
        gap.load_package("guava")
        ans=gap.eval("UpperBoundPlotkin(%s,%s,%s)"%(n,d,q))
        return QQ(ans)
    else:
        t = 1 - 1/q
        if (q==2) and (n == 2*d) and (d%2 == 0):
            return 4*d
        elif (q==2) and (n == 2*d + 1) and (d%2 == 1):
            return 4*d + 4
        elif d > t*n:
            return int(d/( d - t*n))
        elif d < t*n + 1:
            fact = (d-1) / t
            if RR(fact)==RR(int(fact)):
                fact = int(fact) + 1
            return int(d/( d - t * fact)) * q**(n - fact)

def word_problem(words, g, verbose = False):
    r"""
    G and H are abelian, g in G, H is a subgroup of G generated by a
    list (words) of elements of G. If g is in H, return the expression
    for g as a word in the elements of (words).
    
    The 'word problem' for a finite abelian group G boils down to the
    following matrix-vector analog of the Chinese remainder theorem.
    
    Problem: Fix integers `1<n_1\leq n_2\leq ...\leq n_k`
    (indeed, these `n_i` will all be prime powers), fix a
    generating set `g_i=(a_{i1},...,a_{ik})` (with
    `a_{ij}\in \mathrm{Z}/n_j\mathrm{Z}`), for `1\leq i\leq \ell`,
    for the group `G`, and let `d = (d_1,...,d_k)` be
    an element of the direct product
    `\mathrm{Z}/n_1\mathrm{Z} \times ...\times \mathrm{Z}/n_k\mathrm{Z}`. Find, if they
    exist, integers `c_1,...,c_\ell` such that
    `c_1g_1+...+c_\ell g_\ell = d`. In other words, solve
    the equation `cA=d` for `c\in \mathrm{Z}^\ell`, where
    `A` is the matrix whose rows are the `g_i`'s. Of
    course, it suffices to restrict the `c_i`'s to the range
    `0\leq c_i\leq N-1`, where `N` denotes the least
    common multiple of the integers `n_1,...,n_k`.
    
    This function does not solve this directly, as perhaps it should.
    Rather (for both speed and as a model for a similar function valid
    for more general groups), it pushes it over to GAP, which has optimized 
    (non-deterministic) algorithms for the word problem. Essentially, 
    this function is a wrapper for the GAP function 'Factorization'. 
    
    EXAMPLE::
    
        sage: G.<a,b,c> = AbelianGroup(3,[2,3,4]); G
        Multiplicative Abelian Group isomorphic to C2 x C3 x C4
        sage: w = word_problem([a*b,a*c], b*c); w #random 
        [[a*b, 1], [a*c, 1]]  
        sage: prod([x^i for x,i in w]) == b*c
        True
        sage: w = word_problem([a*c,c],a); w #random 
        [[a*c, 1], [c, -1]]   
        sage: prod([x^i for x,i in w]) == a 
        True
        sage: word_problem([a*c,c],a,verbose=True) #random
        a = (a*c)^1*(c)^-1
        [[a*c, 1], [c, -1]]
    
    ::
    
        sage: A.<a,b,c,d,e> = AbelianGroup(5,[4, 5, 5, 7, 8])
        sage: b1 = a^3*b*c*d^2*e^5
        sage: b2 = a^2*b*c^2*d^3*e^3
        sage: b3 = a^7*b^3*c^5*d^4*e^4
        sage: b4 = a^3*b^2*c^2*d^3*e^5
        sage: b5 = a^2*b^4*c^2*d^4*e^5
        sage: w = word_problem([b1,b2,b3,b4,b5],e); w #random
        [[a^3*b*c*d^2*e^5, 1], [a^2*b*c^2*d^3*e^3, 1], [a^3*b^3*d^4*e^4, 3], [a^2*b^4*c^2*d^4*e^5, 1]]
        sage: prod([x^i for x,i in w]) == e 
        True
        sage: word_problem([a,b,c,d,e],e)
        [[e, 1]]
        sage: word_problem([a,b,c,d,e],b)
        [[b, 1]]
    
    .. warning::

       1. Might have unpleasant effect when the word problem
          cannot be solved.
    
       2. Uses permutation groups, so may be slow when group is large.
          The instance method word_problem of the class
          AbelianGroupElement is implemented differently (wrapping
          GAP's 'EpimorphismFromFreeGroup' and
          'PreImagesRepresentative') and may be faster.
       
    """
    from sage.interfaces.all import gap
    G = g.parent()
    invs = G.invariants()
    gap.eval("l:=One(Rationals)")
    s1 = 'A:=AbelianGroup(%s)'%invs
    gap.eval(s1)
    s2 = 'phi:=IsomorphismPermGroup(A)'
    gap.eval(s2)
    s3 = "gens := GeneratorsOfGroup(A)"
    gap.eval(s3)
    L = g.list()
    gap.eval("L1:="+str(L))
    s4 = "L2:=List([l..%s], i->gens[i]^L1[i]);"%len(L) 
    gap.eval(s4)
    gap.eval("g:=Product(L2); gensH:=[]")
    for w in words:
         L = w.list()
         gap.eval("L1:="+str(L))
         s4 = "L2:=List([1..%s], i->gens[i]^L1[i]);"%len(L) 
         gap.eval(s4)
         gap.eval("w:=Product(L2)")
         gap.eval("gensH:=Concatenation(gensH,[w])")
    s5 = 'H:=Group(gensH)'
    gap.eval(s5)
    gap.eval("x:=Factorization(H,g)")
#.........这里部分代码省略.........