def dup_mul(f, g, K):
    """
    Multiply dense polynomials in ``K[x]``.

    Examples
    ========

    >>> from sympy.polys import ring, ZZ
    >>> R, x = ring("x", ZZ)

    >>> R.dup_mul(x - 2, x + 2)
    x**2 - 4

    """
    if f == g:
        return dup_sqr(f, K)

    if not (f and g):
        return []

    df = dup_degree(f)
    dg = dup_degree(g)

    h = []

    for i in xrange(0, df + dg + 1):
        coeff = K.zero

        for j in xrange(max(0, i - dg), min(df, i) + 1):
            coeff += f[j]*g[i - j]

        h.append(coeff)

    return dup_strip(h)

def LU(matlist, K, reverse=0):
    """
    It computes the LU decomposition of a matrix and returns L and U
    matrices.

    Examples
    ========

    >>> from sympy.matrices.densesolve import LU
    >>> from sympy import QQ
    >>> a = [
    ... [QQ(1), QQ(2), QQ(3)],
    ... [QQ(2), QQ(-4), QQ(6)],
    ... [QQ(3), QQ(-9), QQ(-3)]]
    >>> LU(a, QQ)
    ([[1, 0, 0], [2, 1, 0], [3, 15/8, 1]], [[1, 2, 3], [0, -8, 0], [0, 0, -12]])

    See Also
    ========

    upper_triangle
    lower_triangle
    """
    nrow = len(matlist)
    new_matlist1, new_matlist2 = eye(nrow, K), copy.deepcopy(matlist)
    for i in xrange(nrow):
        for j in xrange(i + 1, nrow):
            if new_matlist2[j][i] != 0:
                new_matlist1[j][i] = new_matlist2[j][i] / new_matlist2[i][i]
                rowadd(new_matlist2, j, i, -new_matlist2[j][i] / new_matlist2[i][i], K)
    return new_matlist1, new_matlist2

    def _eval_expand_func(self, **hints):
        n, z = self.args

        if n.is_Integer and n.is_nonnegative:
            if z.is_Add:
                coeff = z.args[0]
                if coeff.is_Integer:
                    e = -(n + 1)
                    if coeff > 0:
                        tail = Add(*[C.Pow(z - i, e) for i in xrange(1, int(coeff) + 1)])
                    else:
                        tail = -Add(*[C.Pow(z + i, e) for i in xrange(0, int(-coeff))])
                    return polygamma(n, z - coeff) + (-1) ** n * C.factorial(n) * tail

            elif z.is_Mul:
                coeff, z = z.as_two_terms()
                if coeff.is_Integer and coeff.is_positive:
                    tail = [polygamma(n, z + C.Rational(i, coeff)) for i in xrange(0, int(coeff))]
                    if n == 0:
                        return Add(*tail) / coeff + log(coeff)
                    else:
                        return Add(*tail) / coeff ** (n + 1)
                z *= coeff

        return polygamma(n, z)

def _identity_matrix(n, domain):
    M = [[domain.zero]*n for _ in xrange(n)]

    for i in xrange(n):
        M[i][i] = domain.one

    return M

def RGS_generalized(m):
    """
    Computes the m + 1 generalized unrestricted growth strings
    and returns them as rows in matrix.

    Examples
    ========

    >>> from sympy.combinatorics.partitions import RGS_generalized
    >>> RGS_generalized(6)
    Matrix([
    [  1,   1,   1,  1,  1, 1, 1],
    [  1,   2,   3,  4,  5, 6, 0],
    [  2,   5,  10, 17, 26, 0, 0],
    [  5,  15,  37, 77,  0, 0, 0],
    [ 15,  52, 151,  0,  0, 0, 0],
    [ 52, 203,   0,  0,  0, 0, 0],
    [203,   0,   0,  0,  0, 0, 0]])
    """
    d = zeros(m + 1)
    for i in xrange(0, m + 1):
        d[0, i] = 1

    for i in xrange(1, m + 1):
        for j in xrange(m):
            if j <= m - i:
                d[i, j] = j * d[i - 1, j] + d[i - 1, j + 1]
            else:
                d[i, j] = 0
    return d

def roots_cyclotomic(f, factor=False):
    """Compute roots of cyclotomic polynomials. """
    L, U = _inv_totient_estimate(f.degree())

    for n in xrange(L, U + 1):
        g = cyclotomic_poly(n, f.gen, polys=True)

        if f == g:
            break
    else:  # pragma: no cover
        raise RuntimeError("failed to find index of a cyclotomic polynomial")

    roots = []

    if not factor:
        # get the indices in the right order so the computed
        # roots will be sorted
        h = n//2
        ks = [i for i in xrange(1, n + 1) if igcd(i, n) == 1]
        ks.sort(key=lambda x: (x, -1) if x <= h else (abs(x - n), 1))
        d = 2*I*pi/n
        for k in reversed(ks):
            roots.append(exp(k*d).expand(complex=True))
    else:
        g = Poly(f, extension=root(-1, n))

        for h, _ in ordered(g.factor_list()[1]):
            roots.append(-h.TC())

    return roots

def eval_levicivita(*args):
    """Evaluate Levi-Civita symbol."""
    from sympy import factorial
    n = len(args)
    return prod(
        prod(args[j] - args[i] for j in xrange(i + 1, n))
        / factorial(i) for i in xrange(n))

def roots_cyclotomic(f, factor=False):
    """Compute roots of cyclotomic polynomials. """
    L, U = _inv_totient_estimate(f.degree())

    for n in xrange(L, U + 1):
        g = cyclotomic_poly(n, f.gen, polys=True)

        if f == g:
            break
    else:  # pragma: no cover
        raise RuntimeError("failed to find index of a cyclotomic polynomial")

    roots = []

    if not factor:
        for k in xrange(1, n + 1):
            if igcd(k, n) == 1:
                roots.append(exp(2*k*S.Pi*I/n).expand(complex=True))
    else:
        g = Poly(f, extension=(-1)**Rational(1, n))

        for h, _ in g.factor_list()[1]:
            roots.append(-h.TC())

    return sorted(roots, key=default_sort_key)

    def _match_common_args(Func, funcs):
        if order != 'none':
            funcs = list(ordered(funcs))
        else:
            funcs = sorted(funcs, key=lambda x: len(x.args))

        func_args = [set(e.args) for e in funcs]
        for i in xrange(len(func_args)):
            for j in xrange(i + 1, len(func_args)):
                com_args = func_args[i].intersection(func_args[j])
                if len(com_args) > 1:
                    com_func = Func(*com_args)

                    # for all sets, replace the common symbols by the function
                    # over them, to allow recursive matches

                    diff_i = func_args[i].difference(com_args)
                    func_args[i] = diff_i | set([com_func])
                    if diff_i:
                        opt_subs[funcs[i]] = Func(Func(*diff_i), com_func,
                                                  evaluate=False)

                    diff_j = func_args[j].difference(com_args)
                    func_args[j] = diff_j | set([com_func])
                    opt_subs[funcs[j]] = Func(Func(*diff_j), com_func,
                                              evaluate=False)

                    for k in xrange(j + 1, len(func_args)):
                        if not com_args.difference(func_args[k]):
                            diff_k = func_args[k].difference(com_args)
                            func_args[k] = diff_k | set([com_func])
                            opt_subs[funcs[k]] = Func(Func(*diff_k), com_func,
                                                      evaluate=False)

    def decrement_part(self, part):
        """Decrements part (a subrange of pstack), if possible, returning
        True iff the part was successfully decremented.

        If you think of the v values in the part as a multi-digit
        integer (least significant digit on the right) this is
        basically decrementing that integer, but with the extra
        constraint that the leftmost digit cannot be decremented to 0.

        Parameters
        ==========

        part
           The part, represented as a list of PartComponent objects,
           which is to be decremented.

        """
        plen = len(part)
        for j in xrange(plen - 1, -1, -1):
            if (j == 0 and part[j].v > 1) or (j > 0 and part[j].v > 0):
                # found val to decrement
                part[j].v -= 1
                # Reset trailing parts back to maximum
                for k in xrange(j + 1, plen):
                    part[k].v = part[k].u
                return True
        return False

    def _initialize_enumeration(self, multiplicities):
        """Allocates and initializes the partition stack.

        This is called from the enumeration/counting routines, so
        there is no need to call it separately."""

        num_components = len(multiplicities)
        # cardinality is the total number of elements, whether or not distinct
        cardinality = sum(multiplicities)

        # pstack is the partition stack, which is segmented by
        # f into parts.
        self.pstack = [PartComponent() for i in
                       xrange(num_components * cardinality + 1)]
        self.f = [0] * (cardinality + 1)

        # Initial state - entire multiset in one part.
        for j in xrange(num_components):
            ps = self.pstack[j]
            ps.c = j
            ps.u = multiplicities[j]
            ps.v = multiplicities[j]

        self.f[0] = 0
        self.f[1] = num_components
        self.lpart = 0

def interpolating_poly(n, x, X='x', Y='y'):
    """Construct Lagrange interpolating polynomial for ``n`` data points. """
    if isinstance(X, str):
        X = symbols("%s:%s" % (X, n))

    if isinstance(Y, str):
        Y = symbols("%s:%s" % (Y, n))

    coeffs = []

    for i in xrange(0, n):
        numer = []
        denom = []

        for j in xrange(0, n):
            if i == j:
                continue

            numer.append(x - X[j])
            denom.append(X[i] - X[j])

        numer = Mul(*numer)
        denom = Mul(*denom)

        coeffs.append(numer/denom)

    return Add(*[ coeff*y for coeff, y in zip(coeffs, Y) ])

def dmp_fateman_poly_F_1(n, K):
    """Fateman's GCD benchmark: trivial GCD """
    u = [K(1), K(0)]

    for i in xrange(0, n):
        u = [dmp_one(i, K), u]

    v = [K(1), K(0), K(0)]

    for i in xrange(0, n):
        v = [dmp_one(i, K), dmp_zero(i), v]

    m = n - 1

    U = dmp_add_term(u, dmp_ground(K(1), m), 0, n, K)
    V = dmp_add_term(u, dmp_ground(K(2), m), 0, n, K)

    f = [[-K(3), K(0)], [], [K(1), K(0), -K(1)]]

    W = dmp_add_term(v, dmp_ground(K(1), m), 0, n, K)
    Y = dmp_raise(f, m, 1, K)

    F = dmp_mul(U, V, n, K)
    G = dmp_mul(W, Y, n, K)

    H = dmp_one(n, K)

    return F, G, H

    def _eval_expand_func(self, **hints):
        from sympy import exp, I, floor, Add, Poly, Dummy, exp_polar, unpolarify
        z, s, a = self.args
        if z == 1:
            return zeta(s, a)
        if s.is_Integer and s <= 0:
            t = Dummy('t')
            p = Poly((t + a)**(-s), t)
            start = 1/(1 - t)
            res = S(0)
            for c in reversed(p.all_coeffs()):
                res += c*start
                start = t*start.diff(t)
            return res.subs(t, z)

        if a.is_Rational:
            # See section 18 of
            #   Kelly B. Roach.  Hypergeometric Function Representations.
            #   In: Proceedings of the 1997 International Symposium on Symbolic and
            #   Algebraic Computation, pages 205-211, New York, 1997. ACM.
            # TODO should something be polarified here?
            add = S(0)
            mul = S(1)
            # First reduce a to the interaval (0, 1]
            if a > 1:
                n = floor(a)
                if n == a:
                    n -= 1
                a -= n
                mul = z**(-n)
                add = Add(*[-z**(k - n)/(a + k)**s for k in xrange(n)])
            elif a <= 0:
                n = floor(-a) + 1
                a += n
                mul = z**n
                add = Add(*[z**(n - 1 - k)/(a - k - 1)**s for k in xrange(n)])

            m, n = S([a.p, a.q])
            zet = exp_polar(2*pi*I/n)
            root = z**(1/n)
            return add + mul*n**(s - 1)*Add(
                *[polylog(s, zet**k*root)._eval_expand_func(**hints)
                  / (unpolarify(zet)**k*root)**m for k in xrange(n)])

        # TODO use minpoly instead of ad-hoc methods when issue 2789 is fixed
        if z.func is exp and (z.args[0]/(pi*I)).is_Rational or z in [-1, I, -I]:
            # TODO reference?
            if z == -1:
                p, q = S([1, 2])
            elif z == I:
                p, q = S([1, 4])
            elif z == -I:
                p, q = S([-1, 4])
            else:
                arg = z.args[0]/(2*pi*I)
                p, q = S([arg.p, arg.q])
            return Add(*[exp(2*pi*I*k*p/q)/q**s*zeta(s, (k + a)/q)
                         for k in xrange(q)])

        return lerchphi(z, s, a)

def dup_from_dict(f, K):
    """
    Create a ``K[x]`` polynomial from a ``dict``.

    Examples
    ========

    >>> from sympy.polys.domains import ZZ
    >>> from sympy.polys.densebasic import dup_from_dict

    >>> dup_from_dict({(0,): ZZ(7), (2,): ZZ(5), (4,): ZZ(1)}, ZZ)
    [1, 0, 5, 0, 7]
    >>> dup_from_dict({}, ZZ)
    []

    """
    if not f:
        return []

    n, h = max(f.keys()), []

    if type(n) is int:
        for k in xrange(n, -1, -1):
            h.append(f.get(k, K.zero))
    else:
        (n,) = n

        for k in xrange(n, -1, -1):
            h.append(f.get((k,), K.zero))

    return dup_strip(h)

def _a(n, j, prec):
    """Compute the inner sum in the HRR formula."""
    if j == 1:
        return fone
    s = fzero
    pi = pi_fixed(prec)
    for h in xrange(1, j):
        if igcd(h, j) != 1:
            continue
        # & with mask to compute fractional part of fixed-point number
        one = 1 << prec
        onemask = one - 1
        half = one >> 1
        g = 0
        if j >= 3:
            for k in xrange(1, j):
                t = h*k*one//j
                if t > 0:
                    frac = t & onemask
                else:
                    frac = -((-t) & onemask)
                g += k*(frac - half)
        g = ((g - 2*h*n*one)*pi//j) >> prec
        s = mpf_add(s, mpf_cos(from_man_exp(g, -prec), prec), prec)
    return s

def _dup_right_decompose(f, s, K):
    """Helper function for :func:`_dup_decompose`."""
    n = len(f) - 1
    lc = dup_LC(f, K)

    f = dup_to_raw_dict(f)
    g = { s: K.one }

    r = n // s

    for i in xrange(1, s):
        coeff = K.zero

        for j in xrange(0, i):
            if not n + j - i in f:
                continue

            if not s - j in g:
                continue

            fc, gc = f[n + j - i], g[s - j]
            coeff += (i - r*j)*fc*gc

        g[s - i] = K.quo(coeff, i*r*lc)

    return dup_from_raw_dict(g, K)

def mdft(n):
    r"""
    Returns an expression of a discrete Fourier transform as a matrix multiplication.
    It is an n X n matrix.

    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/DFT_matrix

    Examples
    ========

    >>> from sympy.physics.matrices import mdft
    >>> mdft(3)
    Matrix([
    [sqrt(3)/3,                sqrt(3)/3,                sqrt(3)/3],
    [sqrt(3)/3, sqrt(3)*exp(-2*I*pi/3)/3, sqrt(3)*exp(-4*I*pi/3)/3],
    [sqrt(3)/3, sqrt(3)*exp(-4*I*pi/3)/3, sqrt(3)*exp(-8*I*pi/3)/3]])
    """
    mat = [[None for x in xrange(n)] for y in xrange(n)]
    base = exp(-2*pi*I/n)
    mat[0] = [1]*n
    for i in range(n):
        mat[i][0] = 1
    for i in xrange(1, n):
        for j in xrange(i, n):
            mat[i][j] = mat[j][i] = base**(i*j)
    return (1/sqrt(n))*Matrix(mat)

    def eval(cls, x, k):
        x = sympify(x)
        k = sympify(k)

        if x is S.NaN:
            return S.NaN
        elif k.is_Integer:
            if k is S.NaN:
                return S.NaN
            elif k is S.Zero:
                return S.One
            else:
                if k.is_positive:
                    if x is S.Infinity:
                        return S.Infinity
                    elif x is S.NegativeInfinity:
                        if k.is_odd:
                            return S.NegativeInfinity
                        else:
                            return S.Infinity
                    else:
                        return reduce(lambda r, i: r*(x - i), xrange(0, int(k)), 1)
                else:
                    if x is S.Infinity:
                        return S.Infinity
                    elif x is S.NegativeInfinity:
                        return S.Infinity
                    else:
                        return 1/reduce(lambda r, i: r*(x + i), xrange(1, abs(int(k)) + 1), 1)

    def _eval_expand_func(self, **hints):
        n = self.args[0]
        m = self.args[1] if len(self.args) == 2 else 1

        if m == S.One:
            if n.is_Add:
                off = n.args[0]
                nnew = n - off
                if off.is_Integer and off.is_positive:
                    result = [S.One/(nnew + i) for i in xrange(off, 0, -1)] + [harmonic(nnew)]
                    return Add(*result)
                elif off.is_Integer and off.is_negative:
                    result = [-S.One/(nnew + i) for i in xrange(0, off, -1)] + [harmonic(nnew)]
                    return Add(*result)

            if n.is_Rational:
                # Expansions for harmonic numbers at general rational arguments (u + p/q)
                # Split n as u + p/q with p < q
                p, q = n.as_numer_denom()
                u = p // q
                p = p - u * q
                if u.is_nonnegative and p.is_positive and q.is_positive and p < q:
                    k = Dummy("k")
                    t1 = q * C.Sum(1 / (q * k + p), (k, 0, u))
                    t2 = 2 * C.Sum(cos((2 * pi * p * k) / S(q)) *
                                   log(sin((pi * k) / S(q))),
                                   (k, 1, floor((q - 1) / S(2))))
                    t3 = (pi / 2) * cot((pi * p) / q) + log(2 * q)
                    return t1 + t2 - t3

        return self