def _solve_as_poly(f, symbol, solveset_solver, invert_func):
    """
    Solve the equation using polynomial techniques if it already is a
    polynomial equation or, with a change of variables, can be made so.
    """
    result = None
    if f.is_polynomial(symbol):

        solns = roots(f, symbol, cubics=True, quartics=True,
                      quintics=True, domain='EX')
        num_roots = sum(solns.values())
        if degree(f, symbol) <= num_roots:
            result = FiniteSet(*solns.keys())
        else:
            poly = Poly(f, symbol)
            solns = poly.all_roots()
            if poly.degree() <= len(solns):
                result = FiniteSet(*solns)
            else:
                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
    else:
        poly = Poly(f)
        if poly is None:
            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
        gens = [g for g in poly.gens if g.has(symbol)]

        if len(gens) == 1:
            poly = Poly(poly, gens[0])
            gen = poly.gen
            deg = poly.degree()
            poly = Poly(poly.as_expr(), poly.gen, composite=True)
            poly_solns = FiniteSet(*roots(poly, cubics=True, quartics=True,
                                          quintics=True).keys())

            if len(poly_solns) < deg:
                result = ConditionSet(symbol, Eq(f, 0), S.Complexes)

            if gen != symbol:
                y = Dummy('y')
                lhs, rhs_s = invert_func(gen, y, symbol)
                if lhs is symbol:
                    result = Union(*[rhs_s.subs(y, s) for s in poly_solns])
                else:
                    result = ConditionSet(symbol, Eq(f, 0), S.Complexes)
        else:
            result = ConditionSet(symbol, Eq(f, 0), S.Complexes)

    if result is not None:
        if isinstance(result, FiniteSet):
            # this is to simplify solutions like -sqrt(-I) to sqrt(2)/2
            # - sqrt(2)*I/2. We are not expanding for solution with free
            # variables because that makes the solution more complicated. For
            # example expand_complex(a) returns re(a) + I*im(a)
            if all([s.free_symbols == set() and not isinstance(s, RootOf)
                    for s in result]):
                s = Dummy('s')
                result = imageset(Lambda(s, expand_complex(s)), result)
        return result
    else:
        return ConditionSet(symbol, Eq(f, 0), S.Complexes)

def _compute_formula(f, x, P, Q, k, m, k_max):
    """Computes the formula for f."""
    from sympy.polys import roots

    sol = []
    for i in range(k_max + 1, k_max + m + 1):
        r = f.diff(x, i).limit(x, 0) / factorial(i)
        if r is S.Zero:
            continue

        kterm = m*k + i
        res = r

        p = P.subs(k, kterm)
        q = Q.subs(k, kterm)
        c1 = p.subs(k, 1/k).leadterm(k)[0]
        c2 = q.subs(k, 1/k).leadterm(k)[0]
        res *= (-c1 / c2)**k

        for r, mul in roots(p, k).items():
            res *= rf(-r, k)**mul
        for r, mul in roots(q, k).items():
            res /= rf(-r, k)**mul

        sol.append((res, kterm))

    return sol

def solve_biquadratic(f, g, opt):
    """Solve a system of two bivariate quadratic polynomial equations. """
    G = groebner([f, g])

    if len(G) == 1 and G[0].is_ground:
        return None

    if len(G) != 2:
        raise SolveFailed

    p, q = G
    x, y = opt.gens

    p = Poly(p, x, expand=False)
    q = q.ltrim(-1)

    p_roots = [ rcollect(expr, y) for expr in roots(p).keys() ]
    q_roots = roots(q).keys()

    solutions = []

    for q_root in q_roots:
        for p_root in p_roots:
            solution = (p_root.subs(y, q_root), q_root)
            solutions.append(solution)

    return sorted(solutions)

    def _solve_reduced_system(system, gens, entry=False):
        """Recursively solves reduced polynomial systems. """
        if len(system) == len(gens) == 1:
            zeros = list(roots(system[0], gens[-1]).keys())
            return [(zero,) for zero in zeros]

        basis = groebner(system, gens, polys=True)

        if len(basis) == 1 and basis[0].is_ground:
            if not entry:
                return []
            else:
                return None

        univariate = list(filter(_is_univariate, basis))

        if len(univariate) == 1:
            f = univariate.pop()
        else:
            raise NotImplementedError(filldedent('''
                only zero-dimensional systems supported
                (finite number of solutions)
                '''))

        gens = f.gens
        gen = gens[-1]

        zeros = list(roots(f.ltrim(gen)).keys())

        if not zeros:
            return []

        if len(basis) == 1:
            return [(zero,) for zero in zeros]

        solutions = []

        for zero in zeros:
            new_system = []
            new_gens = gens[:-1]

            for b in basis[:-1]:
                eq = _subs_root(b, gen, zero)

                if eq is not S.Zero:
                    new_system.append(eq)

            for solution in _solve_reduced_system(new_system, new_gens):
                solutions.append(solution + (zero,))

        if solutions and len(solutions[0]) != len(gens):
            raise NotImplementedError(filldedent('''
                only zero-dimensional systems supported
                (finite number of solutions)
                '''))
        return solutions

def solve_biquadratic(f, g, opt):
    """Solve a system of two bivariate quadratic polynomial equations.

    Examples
    ========

    >>> from sympy.polys import Options, Poly
    >>> from sympy.abc import x, y
    >>> from sympy.solvers.polysys import solve_biquadratic
    >>> NewOption = Options((x, y), {'domain': 'ZZ'})

    >>> a = Poly(y**2 - 4 + x, y, x, domain='ZZ')
    >>> b = Poly(y*2 + 3*x - 7, y, x, domain='ZZ')
    >>> solve_biquadratic(a, b, NewOption)
    [(1/3, 3), (41/27, 11/9)]

    >>> a = Poly(y + x**2 - 3, y, x, domain='ZZ')
    >>> b = Poly(-y + x - 4, y, x, domain='ZZ')
    >>> solve_biquadratic(a, b, NewOption)
    [(7/2 - sqrt(29)/2, -sqrt(29)/2 - 1/2), (sqrt(29)/2 + 7/2, -1/2 + \
      sqrt(29)/2)]
    """
    G = groebner([f, g])

    if len(G) == 1 and G[0].is_ground:
        return None

    if len(G) != 2:
        raise SolveFailed

    x, y = opt.gens
    p, q = G
    if not p.gcd(q).is_ground:
        # not 0-dimensional
        raise SolveFailed

    p = Poly(p, x, expand=False)
    p_roots = [rcollect(expr, y) for expr in roots(p).keys()]

    q = q.ltrim(-1)
    q_roots = list(roots(q).keys())

    solutions = []

    for q_root in q_roots:
        for p_root in p_roots:
            solution = (p_root.subs(y, q_root), q_root)
            solutions.append(solution)

    return sorted(solutions, key=default_sort_key)

    def _eval_product(self, term=None):
        k = self.index
        a = self.lower
        n = self.upper

        if term is None:
            term = self.term

        if not term.has(k):
            return term**(n-a+1)
        elif term.is_polynomial(k):
            poly = term.as_poly(k)

            A = B = Q = S.One
            C_= poly.LC

            all_roots = roots(poly, multiple=True)

            for r in all_roots:
                A *= C.RisingFactorial(a-r, n-a+1)
                Q *= n - r

            if len(all_roots) < poly.degree:
                B = Product(quo(poly, Q.as_poly(k)), (k, a, n))

            return poly.LC**(n-a+1) * A * B
        elif term.is_Add:
            p, q = term.as_numer_denom()

            p = self._eval_product(p)
            q = self._eval_product(q)

            return p / q
        elif term.is_Mul:
            exclude, include = [], []

            for t in term.args:
                p = self._eval_product(t)

                if p is not None:
                    exclude.append(p)
                else:
                    include.append(p)

            if not exclude:
                return None
            else:
                A, B = Mul(*exclude), Mul(*include)
                return A * Product(B, (k, a, n))
        elif term.is_Pow:
            if not term.base.has(k):
                s = sum(term.exp, (k, a, n))

                if not isinstance(s, Sum):
                    return term.base**s
            elif not term.exp.has(k):
                p = self._eval_product(term.base)

                if p is not None:
                    return p**term.exp

    def _eval_product(self, a, n, term):
        from sympy import summation, Sum
        k = self.index

        if not term.has(k):
            return term**(n-a+1)
        elif term.is_polynomial(k):
            poly = term.as_poly(k)

            A = B = Q = S.One
            C_= poly.LC()

            all_roots = roots(poly, multiple=True)

            for r in all_roots:
                A *= C.RisingFactorial(a-r, n-a+1)
                Q *= n - r

            if len(all_roots) < poly.degree():
                B = Product(quo(poly, Q.as_poly(k)), (k, a, n))

            return poly.LC()**(n-a+1) * A * B
        elif term.is_Add:
            p, q = term.as_numer_denom()

            p = self._eval_product(a, n, p)
            q = self._eval_product(a, n, q)

            return p / q
        elif term.is_Mul:
            exclude, include = [], []

            for t in term.args:
                p = self._eval_product(a, n, t)

                if p is not None:
                    exclude.append(p)
                else:
                    include.append(t)

            if not exclude:
                return None
            else:
                A, B = Mul(*exclude), term._new_rawargs(*include)
                return A * Product(B, (k, a, n))
        elif term.is_Pow:
            if not term.base.has(k):
                s = summation(term.exp, (k, a, n))

                if not isinstance(s, Sum):
                    return term.base**s
            elif not term.exp.has(k):
                p = self._eval_product(a, n, term.base)

                if p is not None:
                    return p**term.exp

    def solve_reduced_system(system, gens, entry=False):
        """Recursively solves reduced polynomial systems. """
        basis = groebner(system, gens, polys=True)

        if len(basis) == 1 and basis[0].is_ground:
            if not entry:
                return []
            else:
                return None

        univariate = filter(is_univariate, basis)
        basis = [ b.as_basic() for b in basis ]

        if len(univariate) == 1:
            f = univariate.pop()
        else:
            raise ValueError("only zero-dimensional systems supported")

        gens = f.gens
        f = f.as_basic()

        zeros = roots(f, gens[-1]).keys()

        if not zeros:
            return []

        if len(basis) == 1:
            return [ [zero] for zero in zeros ]

        solutions = []

        for zero in zeros:
            new_system = []
            new_gens = gens[:-1]

            for b in basis[:-1]:
                eq = b.subs(gens[-1], zero).expand()

                if not eq.is_zero:
                    new_system.append(eq)

            for solution in solve_reduced_system(new_system, new_gens):
                solutions.append(solution + [zero])

        return solutions

    def solve_reduced_system(system, entry=False):
        """Recursively solves reduced polynomial systems. """
        basis = poly_groebner(system)

        if len(basis) == 1 and basis[0].is_one:
            if not entry:
                return []
            else:
                return None

        univariate = filter(is_univariate, basis)

        if len(univariate) == 1:
            f = univariate.pop()
        else:
            raise PolynomialError("Not a zero-dimensional system")

        zeros = roots(Poly(f, f.symbols[-1])).keys()

        if not zeros:
            return []

        if len(basis) == 1:
            return [ [zero] for zero in zeros ]

        solutions = []

        for zero in zeros:
            new_system = []

            for poly in basis[:-1]:
                eq = poly.evaluate((poly.symbols[-1], zero))

                if not eq.is_zero:
                    new_system.append(eq)

            for solution in solve_reduced_system(new_system):
                solutions.append(solution + [zero])

        return solutions

#.........这里部分代码省略.........

        # Create a swap dictionary for storing the passed symbols to be solved
        # for, so that they may be swapped back.
        if symbol_swapped:
            swap_dict = zip(symbols, symbols_new)
            f = f.subs(swap_dict)
            symbols = symbols_new

        # Any embedded piecewise functions need to be brought out to the
        # top level so that the appropriate strategy gets selected.
        f = piecewise_fold(f)

        if len(symbols) != 1:
            result = {}
            for s in symbols:
                result[s] = solve(f, s, **flags)
            if flags.get("simplified", True):
                for s, r in result.items():
                    result[s] = map(simplify, r)
            return result

        symbol = symbols[0]
        strategy = guess_solve_strategy(f, symbol)

        if strategy == GS_POLY:
            poly = f.as_poly(symbol)
            if poly is None:
                raise NotImplementedError("Cannot solve equation " + str(f) + " for " + str(symbol))
            # for cubics and quartics, if the flag wasn't set, DON'T do it
            # by default since the results are quite long. Perhaps one could
            # base this decision on a certain crtical length of the roots.
            if poly.degree > 2:
                flags["simplified"] = flags.get("simplified", False)
            result = roots(poly, cubics=True, quartics=True).keys()

        elif strategy == GS_RATIONAL:
            P, Q = f.as_numer_denom()
            # TODO: check for Q != 0
            result = solve(P, symbol, **flags)

        elif strategy == GS_POLY_CV_1:
            args = list(f.args)
            if isinstance(f, Add):
                # we must search for a suitable change of variable
                # collect exponents
                exponents_denom = list()
                for arg in args:
                    if isinstance(arg, Pow):
                        exponents_denom.append(arg.exp.q)
                    elif isinstance(arg, Mul):
                        for mul_arg in arg.args:
                            if isinstance(mul_arg, Pow):
                                exponents_denom.append(mul_arg.exp.q)
                assert len(exponents_denom) > 0
                if len(exponents_denom) == 1:
                    m = exponents_denom[0]
                else:
                    # get the LCM of the denominators
                    m = reduce(ilcm, exponents_denom)
                # x -> y**m.
                # we assume positive for simplification purposes
                t = Dummy("t", positive=True)
                f_ = f.subs(symbol, t ** m)
                if guess_solve_strategy(f_, t) != GS_POLY:
                    raise NotImplementedError("Could not convert to a polynomial equation: %s" % f_)
                cv_sols = solve(f_, t)

def _solve(f, *symbols, **flags):
    """ Return a checked solution for f in terms of one or more of the symbols."""

    if not iterable(f):

        if len(symbols) != 1:
            soln = None
            free = f.free_symbols
            ex = free - set(symbols)
            if len(ex) == 1:
                ex = ex.pop()
                try:
                    # may come back as dict or list (if non-linear)
                    soln = solve_undetermined_coeffs(f, symbols, ex)
                except NotImplementedError:
                    pass
            if not soln is None:
                return soln
            # find first successful solution
            failed = []
            for s in symbols:
                n, d = solve_linear(f, x=[s])
                if n.is_Symbol:
                    soln = {n: cancel(d)}
                    return soln
                failed.append(s)
            for s in failed:
                try:
                    soln = _solve(f, s, **flags)
                    return soln
                except NotImplementedError:
                    pass
            else:
                msg = "No algorithms are implemented to solve equation %s"
                raise NotImplementedError(msg % f)

        symbol = symbols[0]

        # first see if it really depends on symbol and whether there
        # is a linear solution
        f_num, sol = solve_linear(f, x=symbols)
        if not symbol in f_num.free_symbols:
            return []
        elif f_num.is_Symbol:
            return [cancel(sol)]

        strategy = guess_solve_strategy(f, symbol)
        result = False # no solution was obtained

        if strategy == GS_POLY:
            poly = f.as_poly(symbol)
            if poly is None:
                msg = "Cannot solve equation %s for %s" % (f, symbol)
            else:
                # for cubics and quartics, if the flag wasn't set, DON'T do it
                # by default since the results are quite long. Perhaps one could
                # base this decision on a certain critical length of the roots.
                if poly.degree() > 2:
                    flags['simplified'] = flags.get('simplified', False)
                result = roots(poly, cubics=True, quartics=True).keys()

        elif strategy == GS_RATIONAL:
            P, _ = f.as_numer_denom()
            dens = denoms(f, x=symbols)
            try:
                soln = _solve(P, symbol, **flags)
            except NotImplementedError:
                msg = "Cannot solve equation %s for %s" % (P, symbol)
                result = []
            else:
                if dens:
                    # reject any result that makes any denom. affirmatively 0;
                    # if in doubt, keep it
                    result = [s for s in soln if all(not checksol(den, {symbol: s}, **flags) for den in dens)]
                else:
                    result = soln

        elif strategy == GS_POLY_CV_1:
            args = list(f.args)
            if isinstance(f, Pow):
                result = _solve(args[0], symbol, **flags)
            elif isinstance(f, Add):
                # we must search for a suitable change of variables
                # collect exponents
                exponents_denom = list()
                for arg in args:
                    if isinstance(arg, Pow):
                        exponents_denom.append(arg.exp.q)
                    elif isinstance(arg, Mul):
                        for mul_arg in arg.args:
                            if isinstance(mul_arg, Pow):
                                exponents_denom.append(mul_arg.exp.q)
                assert len(exponents_denom) > 0
                if len(exponents_denom) == 1:
                    m = exponents_denom[0]
                else:
                    # get the LCM of the denominators
                    m = reduce(ilcm, exponents_denom)
                # x -> y**m.
                # we assume positive for simplification purposes
#.........这里部分代码省略.........

def solve(f, *symbols, **flags):
    """Solves equations and systems of equations.

       Currently supported are univariate polynomial and transcendental
       equations and systems of linear and polynomial equations.  Input
       is formed as a single expression or an equation,  or an iterable
       container in case of an equation system.  The type of output may
       vary and depends heavily on the input. For more details refer to
       more problem specific functions.

       By default all solutions are simplified to make the output more
       readable. If this is not the expected behavior,  eg. because of
       speed issues, set simplified=False in function arguments.

       To solve equations and systems of equations of other kind, eg.
       recurrence relations of differential equations use rsolve() or
       dsolve() functions respectively.

       >>> from sympy import *
       >>> x,y = symbols('xy')

       Solve a polynomial equation:

       >>> solve(x**4-1, x)
       [1, -1, -I, I]

       Solve a linear system:

       >>> solve((x+5*y-2, -3*x+6*y-15), x, y)
       {x: -3, y: 1}

    """
    if not symbols:
        raise ValueError('no symbols were given')

    if len(symbols) == 1:
        if isinstance(symbols[0], (list, tuple, set)):
            symbols = symbols[0]

    symbols = map(sympify, symbols)
    result = list()

    # Begin code handling for Function and Derivative instances
    # Basic idea:  store all the passed symbols in symbols_passed, check to see
    # if any of them are Function or Derivative types, if so, use a dummy
    # symbol in their place, and set symbol_swapped = True so that other parts
    # of the code can be aware of the swap.  Once all swapping is done, the
    # continue on with regular solving as usual, and swap back at the end of
    # the routine, so that whatever was passed in symbols is what is returned.
    symbols_new = []
    symbol_swapped = False

    if isinstance(symbols, (list, tuple)):
        symbols_passed = symbols[:]
    elif isinstance(symbols, set):
        symbols_passed = list(symbols)

    i = 0
    for s in symbols:
        if s.is_Symbol:
            s_new = s
        elif s.is_Function:
            symbol_swapped = True
            s_new = Symbol('F%d' % i, dummy=True)
        elif s.is_Derivative:
            symbol_swapped = True
            s_new = Symbol('D%d' % i, dummy=True)
        else:
            raise TypeError('not a Symbol or a Function')
        symbols_new.append(s_new)
        i += 1

        if symbol_swapped:
            swap_back_dict = dict(zip(symbols_new, symbols))
    # End code for handling of Function and Derivative instances

    if not isinstance(f, (tuple, list, set)):
        f = sympify(f)

        # Create a swap dictionary for storing the passed symbols to be solved
        # for, so that they may be swapped back.
        if symbol_swapped:
            swap_dict = zip(symbols, symbols_new)
            f = f.subs(swap_dict)
            symbols = symbols_new

        if isinstance(f, Equality):
            f = f.lhs - f.rhs

        if len(symbols) != 1:
            raise NotImplementedError('multivariate equation')

        symbol = symbols[0]

        strategy = guess_solve_strategy(f, symbol)

        if strategy == GS_POLY:
            poly = f.as_poly( symbol )
            assert poly is not None
            result = roots(poly, cubics=True, quartics=True).keys()
#.........这里部分代码省略.........

def rsolve_poly(coeffs, f, n, **hints):
    """Given linear recurrence operator L of order 'k' with polynomial
       coefficients and inhomogeneous equation Ly = f, where 'f' is a
       polynomial, we seek for all polynomial solutions over field K
       of characteristic zero.

       The algorithm performs two basic steps:

           (1) Compute degree N of the general polynomial solution.
           (2) Find all polynomials of degree N or less of Ly = f.

       There are two methods for computing the polynomial solutions.
       If the degree bound is relatively small, i.e. it's smaller than
       or equal to the order of the recurrence, then naive method of
       undetermined coefficients is being used. This gives system
       of algebraic equations with N+1 unknowns.

       In the other case, the algorithm performs transformation of the
       initial equation to an equivalent one, for which the system of
       algebraic equations has only 'r' indeterminates. This method is
       quite sophisticated (in comparison with the naive one) and was
       invented together by Abramov, Bronstein and Petkovsek.

       It is possible to generalize the algorithm implemented here to
       the case of linear q-difference and differential equations.

       Lets say that we would like to compute m-th Bernoulli polynomial
       up to a constant. For this we can use b(n+1) - b(n) == m*n**(m-1)
       recurrence, which has solution b(n) = B_m + C. For example:

       >>> from sympy import Symbol, rsolve_poly
       >>> n = Symbol('n', integer=True)

       >>> rsolve_poly([-1, 1], 4*n**3, n)
       C0 + n**2 - 2*n**3 + n**4

       For more information on implemented algorithms refer to:

       [1] S. A. Abramov, M. Bronstein and M. Petkovsek, On polynomial
           solutions of linear operator equations, in: T. Levelt, ed.,
           Proc. ISSAC '95, ACM Press, New York, 1995, 290-296.

       [2] M. Petkovsek, Hypergeometric solutions of linear recurrences
           with polynomial coefficients, J. Symbolic Computation,
           14 (1992), 243-264.

       [3] M. Petkovsek, H. S. Wilf, D. Zeilberger, A = B, 1996.

    """
    f = sympify(f)

    if not f.is_polynomial(n):
        return None

    homogeneous = f.is_zero

    r = len(coeffs)-1

    coeffs = [ Poly(coeff, n) for coeff in coeffs ]

    polys = [ Poly(0, n) ] * (r+1)
    terms = [ (S.Zero, S.NegativeInfinity) ] *(r+1)

    for i in xrange(0, r+1):
        for j in xrange(i, r+1):
            polys[i] += coeffs[j]*Binomial(j, i)

        if not polys[i].is_zero:
            (exp,), coeff = polys[i].LT()
            terms[i] = (coeff, exp)

    d = b = terms[0][1]

    for i in xrange(1, r+1):
        if terms[i][1] > d:
            d = terms[i][1]

        if terms[i][1] - i > b:
            b = terms[i][1] - i

    d, b = int(d), int(b)

    x = Symbol('x', dummy=True)

    degree_poly = S.Zero

    for i in xrange(0, r+1):
        if terms[i][1] - i == b:
            degree_poly += terms[i][0]*FallingFactorial(x, i)

    nni_roots = roots(degree_poly, x, filter='Z',
        predicate=lambda r: r >= 0).keys()

    if nni_roots:
        N = [max(nni_roots)]
    else:
        N = []

    if homogeneous:
        N += [-b-1]
#.........这里部分代码省略.........

#.........这里部分代码省略.........
            for g, h in similar.iteritems():
                inhomogeneous.append(g+h)
        elif f.is_hypergeometric(n):
            inhomogeneous = [f]
        else:
            return None

        for i, g in enumerate(inhomogeneous):
            coeff, polys = S.One, coeffs[:]
            denoms = [ S.One ] * (r+1)

            s = hypersimp(g, n)

            for j in xrange(1, r+1):
                coeff *= s.subs(n, n+j-1)

                p, q = coeff.as_numer_denom()

                polys[j] *= p
                denoms[j] = q

            for j in xrange(0, r+1):
                polys[j] *= Mul(*(denoms[:j] + denoms[j+1:]))

            R = rsolve_poly(polys, Mul(*denoms), n)

            if not (R is None  or  R is S.Zero):
                inhomogeneous[i] *= R
            else:
                return None

            result = Add(*inhomogeneous)
    else:
        result = S.Zero

    Z = Symbol('Z', dummy=True)

    p, q = coeffs[0], coeffs[r].subs(n, n-r+1)

    p_factors = [ z for z in roots(p, n).iterkeys() ]
    q_factors = [ z for z in roots(q, n).iterkeys() ]

    factors = [ (S.One, S.One) ]

    for p in p_factors:
        for q in q_factors:
            if p.is_integer and q.is_integer and p <= q:
                continue
            else:
                factors += [(n-p, n-q)]

    p = [ (n-p, S.One) for p in p_factors ]
    q = [ (S.One, n-q) for q in q_factors ]

    factors = p + factors + q

    for A, B in factors:
        polys, degrees = [], []
        D = A*B.subs(n, n+r-1)

        for i in xrange(0, r+1):
            a = Mul(*[ A.subs(n, n+j) for j in xrange(0, i) ])
            b = Mul(*[ B.subs(n, n+j) for j in xrange(i, r) ])

            poly = exquo(coeffs[i]*a*b, D, n)
            polys.append(poly.as_poly(n))

            if not poly.is_zero:
                degrees.append(polys[i].degree())

        d, poly = max(degrees), S.Zero

        for i in xrange(0, r+1):
            coeff = polys[i].nth(d)

            if coeff is not S.Zero:
                poly += coeff * Z**i

        for z in roots(poly, Z).iterkeys():
            if not z.is_real or z.is_zero:
                continue

            C = rsolve_poly([ polys[i]*z**i for i in xrange(r+1) ], 0, n)

            if C is not None  and  C is not S.Zero:
                ratio = z * A * C.subs(n, n + 1) / B / C
                K = product(simplify(ratio), (n, 0, n-1))

                if casoratian(kernel+[K], n) != 0:
                    kernel.append(K)

    symbols = [ Symbol('C'+str(i)) for i in xrange(len(kernel)) ]

    for C, ker in zip(symbols, kernel):
        result += C * ker

    if hints.get('symbols', False):
        return (result, symbols)
    else:
        return result

def rsolve_ratio(coeffs, f, n, **hints):
    """Given linear recurrence operator L of order 'k' with polynomial
       coefficients and inhomogeneous equation Ly = f, where 'f' is a
       polynomial, we seek for all rational solutions over field K of
       characteristic zero.

       This procedure accepts only polynomials, however if you are
       interested in solving recurrence with rational coefficients
       then use rsolve() which will pre-process the given equation
       and run this procedure with polynomial arguments.

       The algorithm performs two basic steps:

           (1) Compute polynomial v(n) which can be used as universal
               denominator of any rational solution of equation Ly = f.

           (2) Construct new linear difference equation by substitution
               y(n) = u(n)/v(n) and solve it for u(n) finding all its
               polynomial solutions. Return None if none were found.

       Algorithm implemented here is a revised version of the original
       Abramov's algorithm, developed in 1989. The new approach is much
       simpler to implement and has better overall efficiency. This
       method can be easily adapted to q-difference equations case.

       Besides finding rational solutions alone, this functions is
       an important part of Hyper algorithm were it is used to find
       particular solution of inhomogeneous part of a recurrence.

       For more information on the implemented algorithm refer to:

       [1] S. A. Abramov, Rational solutions of linear difference
           and q-difference equations with polynomial coefficients,
           in: T. Levelt, ed., Proc. ISSAC '95, ACM Press, New York,
           1995, 285-289

    """
    f = sympify(f)

    if not f.is_polynomial(n):
        return None

    coeffs = map(sympify, coeffs)

    r = len(coeffs)-1

    A, B = coeffs[r], coeffs[0]
    A = A.subs(n, n-r).expand()

    h = Symbol('h', dummy=True)

    res = resultant(A, B.subs(n, n+h), n)

    if not res.is_polynomial(h):
        p, q = res.as_numer_denom()
        res = exquo(p, q, h)

    nni_roots = roots(res, h, filter='Z',
        predicate=lambda r: r >= 0).keys()

    if not nni_roots:
        return rsolve_poly(coeffs, f, n, **hints)
    else:
        C, numers = S.One, [S.Zero]*(r+1)

        for i in xrange(int(max(nni_roots)), -1, -1):
            d = gcd(A, B.subs(n, n+i), n)

            A = exquo(A, d, n)
            B = exquo(B, d.subs(n, n-i), n)

            C *= Mul(*[ d.subs(n, n-j) for j in xrange(0, i+1) ])

        denoms = [ C.subs(n, n+i) for i in range(0, r+1) ]

        for i in range(0, r+1):
            g = gcd(coeffs[i], denoms[i], n)

            numers[i] = exquo(coeffs[i], g, n)
            denoms[i] = exquo(denoms[i], g, n)

        for i in xrange(0, r+1):
            numers[i] *= Mul(*(denoms[:i] + denoms[i+1:]))

        result = rsolve_poly(numers, f * Mul(*denoms), n, **hints)

        if result is not None:
            if hints.get('symbols', False):
                return (simplify(result[0] / C), result[1])
            else:
                return simplify(result / C)
        else:
            return None

#.........这里部分代码省略.........

        for i, g in enumerate(inhomogeneous):
            coeff, polys = S.One, coeffs[:]
            denoms = [ S.One ] * (r + 1)

            s = hypersimp(g, n)

            for j in xrange(1, r + 1):
                coeff *= s.subs(n, n + j - 1)

                p, q = coeff.as_numer_denom()

                polys[j] *= p
                denoms[j] = q

            for j in xrange(0, r + 1):
                polys[j] *= Mul(*(denoms[:j] + denoms[j + 1:]))

            R = rsolve_poly(polys, Mul(*denoms), n)

            if not (R is None or R is S.Zero):
                inhomogeneous[i] *= R
            else:
                return None

            result = Add(*inhomogeneous)
    else:
        result = S.Zero

    Z = Dummy('Z')

    p, q = coeffs[0], coeffs[r].subs(n, n - r + 1)

    p_factors = [ z for z in roots(p, n).iterkeys() ]
    q_factors = [ z for z in roots(q, n).iterkeys() ]

    factors = [ (S.One, S.One) ]

    for p in p_factors:
        for q in q_factors:
            if p.is_integer and q.is_integer and p <= q:
                continue
            else:
                factors += [(n - p, n - q)]

    p = [ (n - p, S.One) for p in p_factors ]
    q = [ (S.One, n - q) for q in q_factors ]

    factors = p + factors + q

    for A, B in factors:
        polys, degrees = [], []
        D = A*B.subs(n, n + r - 1)

        for i in xrange(0, r + 1):
            a = Mul(*[ A.subs(n, n + j) for j in xrange(0, i) ])
            b = Mul(*[ B.subs(n, n + j) for j in xrange(i, r) ])

            poly = quo(coeffs[i]*a*b, D, n)
            polys.append(poly.as_poly(n))

            if not poly.is_zero:
                degrees.append(polys[i].degree())

        d, poly = max(degrees), S.Zero

def rsolve_ratio(coeffs, f, n, **hints):
    """
    Given linear recurrence operator `\operatorname{L}` of order `k`
    with polynomial coefficients and inhomogeneous equation
    `\operatorname{L} y = f`, where `f` is a polynomial, we seek
    for all rational solutions over field `K` of characteristic zero.

    This procedure accepts only polynomials, however if you are
    interested in solving recurrence with rational coefficients
    then use ``rsolve`` which will pre-process the given equation
    and run this procedure with polynomial arguments.

    The algorithm performs two basic steps:

        (1) Compute polynomial `v(n)` which can be used as universal
            denominator of any rational solution of equation
            `\operatorname{L} y = f`.

        (2) Construct new linear difference equation by substitution
            `y(n) = u(n)/v(n)` and solve it for `u(n)` finding all its
            polynomial solutions. Return ``None`` if none were found.

    Algorithm implemented here is a revised version of the original
    Abramov's algorithm, developed in 1989. The new approach is much
    simpler to implement and has better overall efficiency. This
    method can be easily adapted to q-difference equations case.

    Besides finding rational solutions alone, this functions is
    an important part of Hyper algorithm were it is used to find
    particular solution of inhomogeneous part of a recurrence.

    Examples
    ========

    >>> from sympy.abc import x
    >>> from sympy.solvers.recurr import rsolve_ratio
    >>> rsolve_ratio([-2*x**3 + x**2 + 2*x - 1, 2*x**3 + x**2 - 6*x,
    ... - 2*x**3 - 11*x**2 - 18*x - 9, 2*x**3 + 13*x**2 + 22*x + 8], 0, x)
    C2*(2*x - 3)/(2*(x**2 - 1))

    References
    ==========

    .. [1] S. A. Abramov, Rational solutions of linear difference
           and q-difference equations with polynomial coefficients,
           in: T. Levelt, ed., Proc. ISSAC '95, ACM Press, New York,
           1995, 285-289

    See Also
    ========

    rsolve_hyper
    """
    f = sympify(f)

    if not f.is_polynomial(n):
        return None

    coeffs = map(sympify, coeffs)

    r = len(coeffs) - 1

    A, B = coeffs[r], coeffs[0]
    A = A.subs(n, n - r).expand()

    h = Dummy('h')

    res = resultant(A, B.subs(n, n + h), n)

    if not res.is_polynomial(h):
        p, q = res.as_numer_denom()
        res = quo(p, q, h)

    nni_roots = roots(res, h, filter='Z',
        predicate=lambda r: r >= 0).keys()

    if not nni_roots:
        return rsolve_poly(coeffs, f, n, **hints)
    else:
        C, numers = S.One, [S.Zero]*(r + 1)

        for i in xrange(int(max(nni_roots)), -1, -1):
            d = gcd(A, B.subs(n, n + i), n)

            A = quo(A, d, n)
            B = quo(B, d.subs(n, n - i), n)

            C *= Mul(*[ d.subs(n, n - j) for j in xrange(0, i + 1) ])

        denoms = [ C.subs(n, n + i) for i in range(0, r + 1) ]

        for i in range(0, r + 1):
            g = gcd(coeffs[i], denoms[i], n)

            numers[i] = quo(coeffs[i], g, n)
            denoms[i] = quo(denoms[i], g, n)

        for i in xrange(0, r + 1):
            numers[i] *= Mul(*(denoms[:i] + denoms[i + 1:]))

#.........这里部分代码省略.........

def solve(f, *symbols, **flags):
    """Solves equations and systems of equations.

       Currently supported are univariate polynomial and transcendental
       equations and systems of linear and polynomial equations.  Input
       is formed as a single expression or an equation,  or an iterable
       container in case of an equation system.  The type of output may
       vary and depends heavily on the input. For more details refer to
       more problem specific functions.

       By default all solutions are simplified to make the output more
       readable. If this is not the expected behavior,  eg. because of
       speed issues, set simplified=False in function arguments.

       To solve equations and systems of equations of other kind, eg.
       recurrence relations of differential equations use rsolve() or
       dsolve() functions respectively.

       >>> from sympy import *
       >>> x,y = symbols('xy')

       Solve a polynomial equation:

       >>> solve(x**4-1, x)
       [1, -1, -I, I]

       Solve a linear system:

       >>> solve((x+5*y-2, -3*x+6*y-15), x, y)
       {x: -3, y: 1}

    """
    if not symbols:
        raise ValueError('no symbols were given')

    if len(symbols) == 1:
        if isinstance(symbols[0], (list, tuple, set)):
            symbols = symbols[0]

    symbols = map(sympify, symbols)

    if any(not s.is_Symbol for s in symbols):
        raise TypeError('not a Symbol')

    if not isinstance(f, (tuple, list, set)):
        f = sympify(f)

        if isinstance(f, Equality):
            f = f.lhs - f.rhs

        if len(symbols) != 1:
            raise NotImplementedError('multivariate equation')

        symbol = symbols[0]

        strategy = guess_solve_stra