def factor_terms(expr, radical=False):
    """Remove common factors from terms in all arguments without
    changing the underlying structure of the expr. No expansion or
    simplification (and no processing of non-commutatives) is performed.

    If radical=True then a radical common to all terms will be factored
    out of any Add sub-expressions of the expr.

    Examples
    ========

    >>> from sympy import factor_terms, Symbol
    >>> from sympy.abc import x, y
    >>> factor_terms(x + x*(2 + 4*y)**3)
    x*(8*(2*y + 1)**3 + 1)
    >>> A = Symbol('A', commutative=False)
    >>> factor_terms(x*A + x*A + x*y*A)
    x*(y*A + 2*A)

    """

    expr = sympify(expr)
    is_iterable = iterable(expr)

    if not isinstance(expr, Basic) or expr.is_Atom:
        if is_iterable:
            return type(expr)([factor_terms(i, radical=radical) for i in expr])
        return expr

    if expr.is_Function or is_iterable or not hasattr(expr, 'args_cnc'):
        args = expr.args
        newargs = tuple([factor_terms(i, radical=radical) for i in args])
        if newargs == args:
            return expr
        return expr.func(*newargs)

    cont, p = expr.as_content_primitive(radical=radical)
    list_args, nc = zip(*[ai.args_cnc() for ai in Add.make_args(p)])
    list_args = list(list_args)
    nc = [((Dummy(), Mul._from_args(i)) if i else None) for i in nc]
    ncreps = dict([i for i in nc if i is not None])
    for i, a in enumerate(list_args):
        if nc[i] is not None:
            a.append(nc[i][0])
        a = Mul._from_args(a) # gcd_terms will fix up ordering
        list_args[i] = gcd_terms(a, isprimitive=True)
        # cancel terms that may not have cancelled
    p = Add._from_args(list_args) # gcd_terms will fix up ordering
    p = gcd_terms(p, isprimitive=True).xreplace(ncreps)
    return _keep_coeff(cont, p)

 def mask(terms):
     """replace nc portions of each term with a unique Dummy symbols
     and return the replacements to restore them"""
     args = [(a, []) if a.is_commutative else a.args_cnc() for a in terms]
     reps = []
     for i, (c, nc) in enumerate(args):
         if nc:
             nc = Mul._from_args(nc)
             d = Dummy()
             reps.append((d, nc))
             c.append(d)
             args[i] = Mul._from_args(c)
         else:
             args[i] = c
     return args, dict(reps)

    def _eval_simplify(self, ratio=1.7, measure=None):
        from sympy.simplify.simplify import factor_sum, sum_combine
        from sympy.core.function import expand
        from sympy.core.mul import Mul

        # split the function into adds
        terms = Add.make_args(expand(self.function))
        s_t = [] # Sum Terms
        o_t = [] # Other Terms

        for term in terms:
            if term.has(Sum):
                # if there is an embedded sum here
                # it is of the form x * (Sum(whatever))
                # hence we make a Mul out of it, and simplify all interior sum terms
                subterms = Mul.make_args(expand(term))
                out_terms = []
                for subterm in subterms:
                    # go through each term
                    if isinstance(subterm, Sum):
                        # if it's a sum, simplify it
                        out_terms.append(subterm._eval_simplify())
                    else:
                        # otherwise, add it as is
                        out_terms.append(subterm)

                # turn it back into a Mul
                s_t.append(Mul(*out_terms))
            else:
                o_t.append(term)

        # next try to combine any interior sums for further simplification
        result = Add(sum_combine(s_t), *o_t)

        return factor_sum(result, limits=self.limits)

def _parse_matrix_expression(expr):
    if isinstance(expr, MatMul):
        args_nonmat = []
        args = []
        contractions = []
        for arg in expr.args:
            if isinstance(arg, MatrixExpr):
                args.append(arg)
            else:
                args_nonmat.append(arg)
        contractions = [(2*i+1, 2*i+2) for i in range(len(args)-1)]
        return Mul.fromiter(args_nonmat)*CodegenArrayContraction(
                CodegenArrayTensorProduct(*[_parse_matrix_expression(arg) for arg in args]),
                *contractions
        )
    elif isinstance(expr, MatAdd):
        return CodegenArrayElementwiseAdd(
                *[_parse_matrix_expression(arg) for arg in expr.args]
        )
    elif isinstance(expr, Transpose):
        return CodegenArrayPermuteDims(
                _parse_matrix_expression(expr.args[0]), [1, 0]
        )
    else:
        return expr

def factor_terms(expr):
    """Remove common factors from terms in all arguments without
    changing the underlying structure of the expr. No expansion or
    simplification (and no processing of non-commutative) is performed.

    **Examples**

    >>> from sympy import factor_terms, Symbol
    >>> from sympy.abc import x, y
    >>> factor_terms(x + x*(2 + 4*y)**3)
    x*(8*(2*y + 1)**3 + 1)
    >>> A = Symbol('A', commutative=False)
    >>> factor_terms(x*A + x*A + x*y*A)
    x*(y*A + 2*A)

    """

    expr = sympify(expr)

    if iterable(expr):
        return type(expr)([factor_terms(i) for i in expr])

    if not isinstance(expr, Basic) or expr.is_Atom:
        return expr

    if expr.is_Function:
        return expr.func(*[factor_terms(i) for i in expr.args])

    cont, p = expr.as_content_primitive()
    list_args, nc = zip(*[ai.args_cnc(clist=True) for ai in Add.make_args(p)])
    list_args = list(list_args)
    nc = [((Dummy(), Mul._from_args(i)) if i else None) for i in nc]
    ncreps = dict([i for i in nc if i is not None])
    for i, a in enumerate(list_args):
        if nc[i] is not None:
           a.append(nc[i][0])
        a = Mul._from_args(a) # gcd_terms will fix up ordering
        list_args[i] = gcd_terms(a, isprimitive=True)
        # cancel terms that may not have cancelled
    p = Add._from_args(list_args) # gcd_terms will fix up ordering
    p = gcd_terms(p, isprimitive=True).subs(ncreps) # exact subs could be used here
    return _keep_coeff(cont, p)

 def from_MatMul(expr):
     args_nonmat = []
     args = []
     contractions = []
     for arg in expr.args:
         if isinstance(arg, MatrixExpr):
             args.append(arg)
         else:
             args_nonmat.append(arg)
     contractions = [(2*i+1, 2*i+2) for i in range(len(args)-1)]
     return Mul.fromiter(args_nonmat)*CodegenArrayContraction(
             CodegenArrayTensorProduct(*args),
             *contractions
         )

 def _sqrt(d):
     # remove squares from square root since both will be represented
     # in the results; a similar thing is happening in roots() but
     # must be duplicated here because not all quadratics are binomials
     co = []
     other = []
     for di in Mul.make_args(d):
         if di.is_Pow and di.exp.is_Integer and di.exp % 2 == 0:
             co.append(Pow(di.base, di.exp//2))
         else:
             other.append(di)
     if co:
         d = Mul(*other)
         co = Mul(*co)
         return co*sqrt(d)
     return sqrt(d)

    def eval(cls, arg):
        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Zero:
                return S.One
            elif arg is S.One:
                return S.Exp1
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.Zero
        elif arg.func is log:
            return arg.args[0]
        elif arg.is_Mul:
            Ioo = S.ImaginaryUnit*S.Infinity
            if arg in [Ioo, -Ioo]:
                return S.NaN

            coeff = arg.coeff(S.Pi*S.ImaginaryUnit)
            if coeff:
                if (2*coeff).is_integer:
                    if coeff.is_even:
                        return S.One
                    elif coeff.is_odd:
                        return S.NegativeOne
                    elif (coeff + S.Half).is_even:
                        return -S.ImaginaryUnit
                    elif (coeff + S.Half).is_odd:
                        return S.ImaginaryUnit

            # Warning: code in risch.py will be very sensitive to changes
            # in this (see DifferentialExtension).

            # look for a single log factor

            coeff, terms = arg.as_coeff_Mul()

            # but it can't be multiplied by oo
            if coeff in [S.NegativeInfinity, S.Infinity]:
                return None

            coeffs, log_term = [coeff], None
            for term in Mul.make_args(terms):
                if term.func is log:
                    if log_term is None:
                        log_term = term.args[0]
                    else:
                        return None
                elif term.is_comparable:
                    coeffs.append(term)
                else:
                    return None

            return log_term**Mul(*coeffs) if log_term else None

        elif arg.is_Add:
            out = []
            add = []
            for a in arg.args:
                if a is S.One:
                    add.append(a)
                    continue
                newa = cls(a)
                if newa.func is cls:
                    add.append(a)
                else:
                    out.append(newa)
            if out:
                return Mul(*out)*cls(Add(*add), evaluate=False)

        elif arg.is_Matrix:
            from sympy import Matrix
            return arg.exp()

def factor_nc(expr):
    """Return the factored form of ``expr`` while handling non-commutative
    expressions.

    **examples**
    >>> from sympy.core.exprtools import factor_nc
    >>> from sympy import Symbol
    >>> from sympy.abc import x
    >>> A = Symbol('A', commutative=False)
    >>> B = Symbol('B', commutative=False)
    >>> factor_nc((x**2 + 2*A*x + A**2).expand())
    (x + A)**2
    >>> factor_nc(((x + A)*(x + B)).expand())
    (x + A)*(x + B)
    """
    from sympy.simplify.simplify import _mexpand
    from sympy.polys import gcd, factor

    expr = sympify(expr)
    if not isinstance(expr, Expr) or not expr.args:
        return expr
    if not expr.is_Add:
        return expr.func(*[factor_nc(a) for a in expr.args])

    expr, rep, nc_symbols = _mask_nc(expr)
    if rep:
        return factor(expr).subs(rep)
    else:
        args = [a.args_cnc() for a in Add.make_args(expr)]
        c = g = l = r = S.One
        hit = False
        # find any commutative gcd term
        for i, a in enumerate(args):
            if i == 0:
                c = Mul._from_args(a[0])
            elif a[0]:
                c = gcd(c, Mul._from_args(a[0]))
            else:
                c = S.One
        if c is not S.One:
            hit = True
            c, g = c.as_coeff_Mul()
            if g is not S.One:
                for i, (cc, _) in enumerate(args):
                    cc = list(Mul.make_args(Mul._from_args(list(cc))/g))
                    args[i][0] = cc
            else:
                for i, (cc, _) in enumerate(args):
                    cc[0] = cc[0]/c
                    args[i][0] = cc
        # find any noncommutative common prefix
        for i, a in enumerate(args):
            if i == 0:
                n = a[1][:]
            else:
                n = common_prefix(n, a[1])
            if not n:
                # is there a power that can be extracted?
                if not args[0][1]:
                    break
                b, e = args[0][1][0].as_base_exp()
                ok = False
                if e.is_Integer:
                    for t in args:
                        if not t[1]:
                            break
                        bt, et = t[1][0].as_base_exp()
                        if et.is_Integer and bt == b:
                            e = min(e, et)
                        else:
                            break
                    else:
                        ok = hit = True
                        l = b**e
                        il = b**-e
                        for i, a in enumerate(args):
                            args[i][1][0] = il*args[i][1][0]
                        break
                if not ok:
                    break
        else:
            hit = True
            lenn = len(n)
            l = Mul(*n)
            for i, a in enumerate(args):
                args[i][1] = args[i][1][lenn:]
        # find any noncommutative common suffix
        for i, a in enumerate(args):
            if i == 0:
                n = a[1][:]
            else:
                n = common_suffix(n, a[1])
            if not n:
                # is there a power that can be extracted?
                if not args[0][1]:
                    break
                b, e = args[0][1][-1].as_base_exp()
                ok = False
                if e.is_Integer:
                    for t in args:
#.........这里部分代码省略.........

    def __init__(self, factors=None):  # Factors
        """Initialize Factors from dict or expr.

        Examples
        ========

        >>> from sympy.core.exprtools import Factors
        >>> from sympy.abc import x
        >>> from sympy import I
        >>> e = 2*x**3
        >>> Factors(e)
        Factors({2: 1, x: 3})
        >>> Factors(e.as_powers_dict())
        Factors({2: 1, x: 3})
        >>> f = _
        >>> f.factors  # underlying dictionary
        {2: 1, x: 3}
        >>> f.gens  # base of each factor
        frozenset([2, x])
        >>> Factors(0)
        Factors({0: 1})
        >>> Factors(I)
        Factors({I: 1})

        Notes
        =====

        Although a dictionary can be passed, only minimal checking is
        performed: powers of -1 and I are made canonical.

        """
        if isinstance(factors, (SYMPY_INTS, float)):
            factors = S(factors)

        if isinstance(factors, Factors):
            factors = factors.factors.copy()
        elif factors is None or factors is S.One:
            factors = {}
        elif factors is S.Zero or factors == 0:
            factors = {S.Zero: S.One}
        elif isinstance(factors, Number):
            n = factors
            factors = {}
            if n < 0:
                factors[S.NegativeOne] = S.One
                n = -n
            if n is not S.One:
                if n.is_Float or n.is_Integer or n is S.Infinity:
                    factors[n] = S.One
                elif n.is_Rational:
                    # since we're processing Numbers, the denominator is
                    # stored with a negative exponent; all other factors
                    # are left .
                    if n.p != 1:
                        factors[Integer(n.p)] = S.One
                    factors[Integer(n.q)] = S.NegativeOne
                else:
                    raise ValueError('Expected Float|Rational|Integer, not %s' % n)
        elif isinstance(factors, Basic) and not factors.args:
            factors = {factors: S.One}
        elif isinstance(factors, Expr):
            c, nc = factors.args_cnc()
            i = c.count(I)
            for _ in range(i):
                c.remove(I)
            factors = dict(Mul._from_args(c).as_powers_dict())
            if i:
                factors[I] = S.One*i
            if nc:
                factors[Mul(*nc, evaluate=False)] = S.One
        else:
            factors = factors.copy()  # /!\ should be dict-like

            # tidy up -/+1 and I exponents if Rational

            handle = []
            for k in factors:
                if k is I or k in (-1, 1):
                    handle.append(k)
            if handle:
                i1 = S.One
                for k in handle:
                    if not _isnumber(factors[k]):
                        continue
                    i1 *= k**factors.pop(k)
                if i1 is not S.One:
                    for a in i1.args if i1.is_Mul else [i1]:  # at worst, -1.0*I*(-1)**e
                        if a is S.NegativeOne:
                            factors[a] = S.One
                        elif a is I:
                            factors[I] = S.One
                        elif a.is_Pow:
                            if S.NegativeOne not in factors:
                                factors[S.NegativeOne] = S.Zero
                            factors[S.NegativeOne] += a.exp
                        elif a == 1:
                            factors[a] = S.One
                        elif a == -1:
                            factors[-a] = S.One
                            factors[S.NegativeOne] = S.One
#.........这里部分代码省略.........

def rsolve(f, y, init=None):
    """
    Solve univariate recurrence with rational coefficients.

    Given `k`-th order linear recurrence `\operatorname{L} y = f`,
    or equivalently:

    .. math:: a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) +
              \dots + a_{0}(n) y(n) = f(n)

    where `a_{i}(n)`, for `i=0, \dots, k`, are polynomials or rational
    functions in `n`, and `f` is a hypergeometric function or a sum
    of a fixed number of pairwise dissimilar hypergeometric terms in
    `n`, finds all solutions or returns ``None``, if none were found.

    Initial conditions can be given as a dictionary in two forms:

        (1) ``{   n_0  : v_0,   n_1  : v_1, ...,   n_m  : v_m }``
        (2) ``{ y(n_0) : v_0, y(n_1) : v_1, ..., y(n_m) : v_m }``

    or as a list ``L`` of values:

        ``L = [ v_0, v_1, ..., v_m ]``

    where ``L[i] = v_i``, for `i=0, \dots, m`, maps to `y(n_i)`.

    Examples
    ========

    Lets consider the following recurrence:

    .. math:: (n - 1) y(n + 2) - (n^2 + 3 n - 2) y(n + 1) +
              2 n (n + 1) y(n) = 0

    >>> from sympy import Function, rsolve
    >>> from sympy.abc import n
    >>> y = Function('y')

    >>> f = (n - 1)*y(n + 2) - (n**2 + 3*n - 2)*y(n + 1) + 2*n*(n + 1)*y(n)

    >>> rsolve(f, y(n))
    2**n*C0 + C1*factorial(n)

    >>> rsolve(f, y(n), { y(0):0, y(1):3 })
    3*2**n - 3*factorial(n)

    See Also
    ========

    rsolve_poly, rsolve_ratio, rsolve_hyper

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

    n = y.args[0]
    k = Wild('k', exclude=(n,))

    # Preprocess user input to allow things like
    # y(n) + a*(y(n + 1) + y(n - 1))/2
    f = f.expand().collect(y.func(Wild('m', integer=True)))

    h_part = defaultdict(lambda: S.Zero)
    i_part = S.Zero
    for g in Add.make_args(f):
        coeff = S.One
        kspec = None
        for h in Mul.make_args(g):
            if h.is_Function:
                if h.func == y.func:
                    result = h.args[0].match(n + k)

                    if result is not None:
                        kspec = int(result[k])
                    else:
                        raise ValueError(
                            "'%s(%s+k)' expected, got '%s'" % (y.func, n, h))
                else:
                    raise ValueError(
                        "'%s' expected, got '%s'" % (y.func, h.func))
            else:
                coeff *= h

        if kspec is not None:
            h_part[kspec] += coeff
        else:
            i_part += coeff

    for k, coeff in h_part.iteritems():
        h_part[k] = simplify(coeff)

    common = S.One

    for coeff in h_part.itervalues():
        if coeff.is_rational_function(n):
            if not coeff.is_polynomial(n):
                common = lcm(common, coeff.as_numer_denom()[1], n)
        else:
            raise ValueError(
                "Polynomial or rational function expected, got '%s'" % coeff)
#.........这里部分代码省略.........

    def eval(cls, p, q):
        from sympy.core.add import Add
        from sympy.core.mul import Mul
        from sympy.core.singleton import S
        from sympy.core.exprtools import gcd_terms
        from sympy.polys.polytools import gcd

        def doit(p, q):
            """Try to return p % q if both are numbers or +/-p is known
            to be less than or equal q.
            """

            if q == S.Zero:
                raise ZeroDivisionError("Modulo by zero")
            if p.is_infinite or q.is_infinite or p is nan or q is nan:
                return nan
            if p == S.Zero or p == q or p == -q or (p.is_integer and q == 1):
                return S.Zero

            if q.is_Number:
                if p.is_Number:
                    return p%q
                if q == 2:
                    if p.is_even:
                        return S.Zero
                    elif p.is_odd:
                        return S.One

            if hasattr(p, '_eval_Mod'):
                rv = getattr(p, '_eval_Mod')(q)
                if rv is not None:
                    return rv

            # by ratio
            r = p/q
            try:
                d = int(r)
            except TypeError:
                pass
            else:
                if isinstance(d, integer_types):
                    rv = p - d*q
                    if (rv*q < 0) == True:
                        rv += q
                    return rv

            # by difference
            # -2|q| < p < 2|q|
            d = abs(p)
            for _ in range(2):
                d -= abs(q)
                if d.is_negative:
                    if q.is_positive:
                        if p.is_positive:
                            return d + q
                        elif p.is_negative:
                            return -d
                    elif q.is_negative:
                        if p.is_positive:
                            return d
                        elif p.is_negative:
                            return -d + q
                    break

        rv = doit(p, q)
        if rv is not None:
            return rv

        # denest
        if isinstance(p, cls):
            qinner = p.args[1]
            if qinner % q == 0:
                return cls(p.args[0], q)
            elif (qinner*(q - qinner)).is_nonnegative:
                # |qinner| < |q| and have same sign
                return p
        elif isinstance(-p, cls):
            qinner = (-p).args[1]
            if qinner % q == 0:
                return cls(-(-p).args[0], q)
            elif (qinner*(q + qinner)).is_nonpositive:
                # |qinner| < |q| and have different sign
                return p
        elif isinstance(p, Add):
            # separating into modulus and non modulus
            both_l = non_mod_l, mod_l = [], []
            for arg in p.args:
                both_l[isinstance(arg, cls)].append(arg)
            # if q same for all
            if mod_l and all(inner.args[1] == q for inner in mod_l):
                net = Add(*non_mod_l) + Add(*[i.args[0] for i in mod_l])
                return cls(net, q)

        elif isinstance(p, Mul):
            # separating into modulus and non modulus
            both_l = non_mod_l, mod_l = [], []
            for arg in p.args:
                both_l[isinstance(arg, cls)].append(arg)

            if mod_l and all(inner.args[1] == q for inner in mod_l):
#.........这里部分代码省略.........

#.........这里部分代码省略.........
            from sympy.polys.polytools import real_roots
            from sympy.polys.polyroots import roots
            from sympy.polys.polyerrors import PolynomialError
            x = free.pop()
            x0 = _monotonic_sign(x)
            if x0 == _eps or x0 == -_eps:
                x0 = S.Zero
            if x0 is not None:
                d = self.diff(x)
                if d.is_number:
                    roots = []
                else:
                    try:
                        roots = real_roots(d)
                    except (PolynomialError, NotImplementedError):
                        roots = [r for r in roots(d, x) if r.is_real]
                y = self.subs(x, x0)
                if x.is_nonnegative and all(r <= x0 for r in roots):
                    if y.is_nonnegative and d.is_positive:
                        if y:
                            return y if y.is_positive else Dummy('pos', positive=True)
                        else:
                            return Dummy('nneg', nonnegative=True)
                    if y.is_nonpositive and d.is_negative:
                        if y:
                            return y if y.is_negative else Dummy('neg', negative=True)
                        else:
                            return Dummy('npos', nonpositive=True)
                elif x.is_nonpositive and all(r >= x0 for r in roots):
                    if y.is_nonnegative and d.is_negative:
                        if y:
                            return Dummy('pos', positive=True)
                        else:
                            return Dummy('nneg', nonnegative=True)
                    if y.is_nonpositive and d.is_positive:
                        if y:
                            return Dummy('neg', negative=True)
                        else:
                            return Dummy('npos', nonpositive=True)
        else:
            n, d = self.as_numer_denom()
            den = None
            if n.is_number:
                den = _monotonic_sign(d)
            elif not d.is_number:
                if _monotonic_sign(n) is not None:
                    den = _monotonic_sign(d)
            if den is not None and (den.is_positive or den.is_negative):
                v = n*den
                if v.is_positive:
                    return Dummy('pos', positive=True)
                elif v.is_nonnegative:
                    return Dummy('nneg', nonnegative=True)
                elif v.is_negative:
                    return Dummy('neg', negative=True)
                elif v.is_nonpositive:
                    return Dummy('npos', nonpositive=True)
        return None

    # multivariate
    c, a = self.as_coeff_Add()
    v = None
    if not a.is_polynomial():
        # F/A or A/F where A is a number and F is a signed, rational monomial
        n, d = a.as_numer_denom()
        if not (n.is_number or d.is_number):
            return
        if (
                a.is_Mul or a.is_Pow) and \
                a.is_rational and \
                all(p.exp.is_Integer for p in a.atoms(Pow) if p.is_Pow) and \
                (a.is_positive or a.is_negative):
            v = S(1)
            for ai in Mul.make_args(a):
                if ai.is_number:
                    v *= ai
                    continue
                reps = {}
                for x in ai.free_symbols:
                    reps[x] = _monotonic_sign(x)
                    if reps[x] is None:
                        return
                v *= ai.subs(reps)
    elif c:
        # signed linear expression
        if not any(p for p in a.atoms(Pow) if not p.is_number) and (a.is_nonpositive or a.is_nonnegative):
            free = list(a.free_symbols)
            p = {}
            for i in free:
                v = _monotonic_sign(i)
                if v is None:
                    return
                p[i] = v or (_eps if i.is_nonnegative else -_eps)
            v = a.xreplace(p)
    if v is not None:
        rv = v + c
        if v.is_nonnegative and rv.is_positive:
            return rv.subs(_eps, 0)
        if v.is_nonpositive and rv.is_negative:
            return rv.subs(_eps, 0)

    def eval(cls, arg):
        from sympy.assumptions import ask, Q
        from sympy.calculus import AccumBounds
        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Zero:
                return S.One
            elif arg is S.One:
                return S.Exp1
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.Zero
        elif isinstance(arg, log):
            return arg.args[0]
        elif isinstance(arg, AccumBounds):
            return AccumBounds(exp(arg.min), exp(arg.max))
        elif arg.is_Mul:
            if arg.is_number or arg.is_Symbol:
                coeff = arg.coeff(S.Pi*S.ImaginaryUnit)
                if coeff:
                    if ask(Q.integer(2*coeff)):
                        if ask(Q.even(coeff)):
                            return S.One
                        elif ask(Q.odd(coeff)):
                            return S.NegativeOne
                        elif ask(Q.even(coeff + S.Half)):
                            return -S.ImaginaryUnit
                        elif ask(Q.odd(coeff + S.Half)):
                            return S.ImaginaryUnit

            # Warning: code in risch.py will be very sensitive to changes
            # in this (see DifferentialExtension).

            # look for a single log factor

            coeff, terms = arg.as_coeff_Mul()

            # but it can't be multiplied by oo
            if coeff in [S.NegativeInfinity, S.Infinity]:
                return None

            coeffs, log_term = [coeff], None
            for term in Mul.make_args(terms):
                if isinstance(term, log):
                    if log_term is None:
                        log_term = term.args[0]
                    else:
                        return None
                elif term.is_comparable:
                    coeffs.append(term)
                else:
                    return None

            return log_term**Mul(*coeffs) if log_term else None

        elif arg.is_Add:
            out = []
            add = []
            for a in arg.args:
                if a is S.One:
                    add.append(a)
                    continue
                newa = cls(a)
                if isinstance(newa, cls):
                    add.append(a)
                else:
                    out.append(newa)
            if out:
                return Mul(*out)*cls(Add(*add), evaluate=False)

        elif arg.is_Matrix:
            return arg.exp()

def rsolve(f, y, init=None):
    """
    Solve univariate recurrence with rational coefficients.

    Given k-th order linear recurrence Ly = f, or equivalently:

     a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) + ... + a_{0}(n) y(n) = f

    where a_{i}(n), for i=0..k, are polynomials or rational functions
    in n, and f is a hypergeometric function or a sum of a fixed number
    of pairwise dissimilar hypergeometric terms in n, finds all solutions
    or returns None, if none were found.

    Initial conditions can be given as a dictionary in two forms:

      [1] {   n_0  : v_0,   n_1  : v_1, ...,   n_m  : v_m }
      [2] { y(n_0) : v_0, y(n_1) : v_1, ..., y(n_m) : v_m }

    or as a list L of values:

      L = [ v_0, v_1, ..., v_m ]

    where L[i] = v_i, for i=0..m, maps to y(n_i).

    As an example lets consider the following recurrence:

     (n - 1) y(n + 2) - (n**2 + 3 n - 2) y(n + 1) + 2 n (n + 1) y(n) == 0

    >>> from sympy import Function, rsolve
    >>> from sympy.abc import n
    >>> y = Function('y')

    >>> f = (n-1)*y(n+2) - (n**2+3*n-2)*y(n+1) + 2*n*(n+1)*y(n)

    >>> rsolve(f, y(n))
    2**n*C0 + C1*n!

    >>> rsolve(f, y(n), { y(0):0, y(1):3 })
    3*2**n - 3*n!

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

    n = y.args[0]
    k = Wild('k', exclude=(n,))

    h_part = defaultdict(lambda: S.Zero)
    i_part = S.Zero
    for g in Add.make_args(f):
        coeff = S.One
        kspec = None
        for h in Mul.make_args(g):
            if h.is_Function:
                if h.func == y.func:
                    result = h.args[0].match(n + k)

                    if result is not None:
                        kspec = int(result[k])
                    else:
                        raise ValueError(
                            "'%s(%s+k)' expected, got '%s'" % (y.func, n, h))
                else:
                    raise ValueError(
                        "'%s' expected, got '%s'" % (y.func, h.func))
            else:
                coeff *= h

        if kspec is not None:
            h_part[kspec] += coeff
        else:
            i_part += coeff

    for k, coeff in h_part.iteritems():
        h_part[k] = simplify(coeff)

    common = S.One

    for coeff in h_part.itervalues():
        if coeff.is_rational_function(n):
            if not coeff.is_polynomial(n):
                common = lcm(common, coeff.as_numer_denom()[1], n)
        else:
            raise ValueError(
                "Polynomial or rational function expected, got '%s'" % coeff)

    i_numer, i_denom = i_part.as_numer_denom()

    if i_denom.is_polynomial(n):
        common = lcm(common, i_denom, n)

    if common is not S.One:
        for k, coeff in h_part.iteritems():
            numer, denom = coeff.as_numer_denom()
            h_part[k] = numer*quo(common, denom, n)

        i_part = i_numer*quo(common, i_denom, n)

    K_min = min(h_part.keys())

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

    def eval(cls, p, q):
        from sympy.core.add import Add
        from sympy.core.mul import Mul
        from sympy.core.singleton import S
        from sympy.core.exprtools import gcd_terms
        from sympy.polys.polytools import gcd

        def doit(p, q):
            """Try to return p % q if both are numbers or +/-p is known
            to be less than q.
            """

            if p == q or p == -q or p.is_Pow and p.exp.is_Integer and p.base == q:
                return S.Zero

            if p.is_Number and q.is_Number:
                return (p % q)

            # by ratio
            r = p/q
            try:
                d = int(r)
            except TypeError:
                pass
            else:
                if type(d) is int:
                    rv = p - d*q
                    if (rv*q < 0) is True:
                        rv += q
                    return rv

            # by differencec
            d = p - q
            if d.is_negative:
                if q.is_negative:
                    return d
                elif q.is_positive:
                    return p

        rv = doit(p, q)
        if rv is not None:
            return rv

        # denest
        if p.func is cls:
            # easy
            qinner = p.args[1]
            if qinner == q:
                return p
            # XXX other possibilities?

        # extract gcd; any further simplification should be done by the user
        G = gcd(p, q)
        if G is not S.One:
            p, q = [
                gcd_terms(i/G, clear=False, fraction=False) for i in (p, q)]
        pwas, qwas = p, q

        # simplify terms
        # (x + y + 2) % x -> Mod(y + 2, x)
        if p.is_Add:
            args = []
            for i in p.args:
                a = cls(i, q)
                if a.count(cls) > i.count(cls):
                    args.append(i)
                else:
                    args.append(a)
            if args != list(p.args):
                p = Add(*args)

        else:
            # handle coefficients if they are not Rational
            # since those are not handled by factor_terms
            # e.g. Mod(.6*x, .3*y) -> 0.3*Mod(2*x, y)
            cp, p = p.as_coeff_Mul()
            cq, q = q.as_coeff_Mul()
            ok = False
            if not cp.is_Rational or not cq.is_Rational:
                r = cp % cq
                if r == 0:
                    G *= cq
                    p *= int(cp/cq)
                    ok = True
            if not ok:
                p = cp*p
                q = cq*q

        # simple -1 extraction
        if p.could_extract_minus_sign() and q.could_extract_minus_sign():
            G, p, q = [-i for i in (G, p, q)]

        # check again to see if p and q can now be handled as numbers
        rv = doit(p, q)
        if rv is not None:
            return rv*G

        # put 1.0 from G on inside
        if G.is_Float and G == 1:
            p *= G
#.........这里部分代码省略.........

def roots(f, *gens, **flags):
    """
    Computes symbolic roots of a univariate polynomial.

    Given a univariate polynomial f with symbolic coefficients (or
    a list of the polynomial's coefficients), returns a dictionary
    with its roots and their multiplicities.

    Only roots expressible via radicals will be returned.  To get
    a complete set of roots use RootOf class or numerical methods
    instead. By default cubic and quartic formulas are used in
    the algorithm. To disable them because of unreadable output
    set ``cubics=False`` or ``quartics=False`` respectively. If cubic
    roots are real but are expressed in terms of complex numbers
    (casus irreducibilis [1]) the ``trig`` flag can be set to True to
    have the solutions returned in terms of cosine and inverse cosine
    functions.

    To get roots from a specific domain set the ``filter`` flag with
    one of the following specifiers: Z, Q, R, I, C. By default all
    roots are returned (this is equivalent to setting ``filter='C'``).

    By default a dictionary is returned giving a compact result in
    case of multiple roots.  However to get a list containing all
    those roots set the ``multiple`` flag to True; the list will
    have identical roots appearing next to each other in the result.
    (For a given Poly, the all_roots method will give the roots in
    sorted numerical order.)

    Examples
    ========

    >>> from sympy import Poly, roots
    >>> from sympy.abc import x, y

    >>> roots(x**2 - 1, x)
    {-1: 1, 1: 1}

    >>> p = Poly(x**2-1, x)
    >>> roots(p)
    {-1: 1, 1: 1}

    >>> p = Poly(x**2-y, x, y)

    >>> roots(Poly(p, x))
    {-sqrt(y): 1, sqrt(y): 1}

    >>> roots(x**2 - y, x)
    {-sqrt(y): 1, sqrt(y): 1}

    >>> roots([1, 0, -1])
    {-1: 1, 1: 1}


    References
    ==========

    .. [1] https://en.wikipedia.org/wiki/Cubic_function#Trigonometric_.28and_hyperbolic.29_method

    """
    from sympy.polys.polytools import to_rational_coeffs
    flags = dict(flags)

    auto = flags.pop('auto', True)
    cubics = flags.pop('cubics', True)
    trig = flags.pop('trig', False)
    quartics = flags.pop('quartics', True)
    quintics = flags.pop('quintics', False)
    multiple = flags.pop('multiple', False)
    filter = flags.pop('filter', None)
    predicate = flags.pop('predicate', None)

    if isinstance(f, list):
        if gens:
            raise ValueError('redundant generators given')

        x = Dummy('x')

        poly, i = {}, len(f) - 1

        for coeff in f:
            poly[i], i = sympify(coeff), i - 1

        f = Poly(poly, x, field=True)
    else:
        try:
            f = Poly(f, *gens, **flags)
            if f.length == 2 and f.degree() != 1:
                # check for foo**n factors in the constant
                n = f.degree()
                npow_bases = []
                others = []
                expr = f.as_expr()
                con = expr.as_independent(*gens)[0]
                for p in Mul.make_args(con):
                    if p.is_Pow and not p.exp % n:
                        npow_bases.append(p.base**(p.exp/n))
                    else:
                        others.append(p)
                    if npow_bases:
#.........这里部分代码省略.........