def _mul(self, other):
        """See docstring of SeqBase._mul"""
        if isinstance(other, SeqPer):
            per1, lper1 = self.periodical, self.period
            per2, lper2 = other.periodical, other.period

            per_length = lcm(lper1, lper2)

            new_per = []
            for x in range(per_length):
                ele1 = per1[x % lper1]
                ele2 = per2[x % lper2]
                new_per.append(ele1 * ele2)

            start, stop = self._intersect_interval(other)
            return SeqPer(new_per, (self.variables[0], start, stop))

def _rsolve_hypergeometric(f, x, P, Q, k, m):
    """Recursive wrapper to rsolve_hypergeometric.

    Returns a Tuple of (formula, series independent terms,
    maximum power of x in independent terms) if successful
    otherwise ``None``.

    See :func:`rsolve_hypergeometric` for details.
    """
    from sympy.polys import lcm, roots
    from sympy.integrals import integrate

    # tranformation - c
    proots, qroots = roots(P, k), roots(Q, k)
    all_roots = dict(proots)
    all_roots.update(qroots)
    scale = lcm([r.as_numer_denom()[1] for r, t in all_roots.items()
                 if r.is_rational])
    f, P, Q, m = _transformation_c(f, x, P, Q, k, m, scale)

    # transformation - a
    qroots = roots(Q, k)
    if qroots:
        k_min = Min(*qroots.keys())
    else:
        k_min = S.Zero
    shift = k_min + m
    f, P, Q, m = _transformation_a(f, x, P, Q, k, m, shift)

    l = (x*f).limit(x, 0)
    if not isinstance(l, Limit) and l != 0:  # Ideally should only be l != 0
        return None

    qroots = roots(Q, k)
    if qroots:
        k_max = Max(*qroots.keys())
    else:
        k_max = S.Zero

    ind, mp = S.Zero, -oo
    for i in range(k_max + m + 1):
        r = f.diff(x, i).limit(x, 0) / factorial(i)
        if r.is_finite is False:
            old_f = f
            f, P, Q, m = _transformation_a(f, x, P, Q, k, m, i)
            f, P, Q, m = _transformation_e(f, x, P, Q, k, m)
            sol, ind, mp = _rsolve_hypergeometric(f, x, P, Q, k, m)
            sol = _apply_integrate(sol, x, k)
            sol = _apply_shift(sol, i)
            ind = integrate(ind, x)
            ind += (old_f - ind).limit(x, 0)  # constant of integration
            mp += 1
            return sol, ind, mp
        elif r:
            ind += r*x**(i + shift)
            pow_x = Rational((i + shift), scale)
            if pow_x > mp:
                mp = pow_x  # maximum power of x
    ind = ind.subs(x, x**(1/scale))

    sol = _compute_formula(f, x, P, Q, k, m, k_max)
    sol = _apply_shift(sol, shift)
    sol = _apply_scale(sol, scale)

    return sol, ind, mp

#.........这里部分代码省略.........
        terms |= components(cancel(g.diff(x)), x)

    # TODO: caching is significant factor for why permutations work at all. Change this.
    V = _symbols('x', len(terms))

    mapping = dict(list(zip(terms, V)))

    rev_mapping = {}

    if unnecessary_permutations is None:
        unnecessary_permutations = []
    for k, v in mapping.items():
        rev_mapping[v] = k

    if mappings is None:
        # Pre-sort mapping in order of largest to smallest expressions (last is always x).
        def _sort_key(arg):
            return default_sort_key(arg[0].as_independent(x)[1])
        #optimizing the number of permutations of mappping
        unnecessary_permutations = [(x, mapping[x])]
        del mapping[x]
        mapping = sorted(list(mapping.items()), key=_sort_key, reverse=True)
        mappings = permutations(mapping)

    def _substitute(expr):
        return expr.subs(mapping)

    for mapping in mappings:
        mapping = list(mapping)
        mapping = mapping + unnecessary_permutations
        diffs = [ _substitute(cancel(g.diff(x))) for g in terms ]
        denoms = [ g.as_numer_denom()[1] for g in diffs ]
        if all(h.is_polynomial(*V) for h in denoms) and _substitute(f).is_rational_function(*V):
            denom = reduce(lambda p, q: lcm(p, q, *V), denoms)
            break
    else:
        if not rewrite:
            result = heurisch(f, x, rewrite=True, hints=hints, unnecessary_permutations=unnecessary_permutations)

            if result is not None:
                return indep*result
        return None

    numers = [ cancel(denom*g) for g in diffs ]
    def _derivation(h):
        return Add(*[ d * h.diff(v) for d, v in zip(numers, V) ])

    def _deflation(p):
        for y in V:
            if not p.has(y):
                continue

            if _derivation(p) is not S.Zero:
                c, q = p.as_poly(y).primitive()
                return _deflation(c)*gcd(q, q.diff(y)).as_expr()
        else:
            return p

    def _splitter(p):
        for y in V:
            if not p.has(y):
                continue

            if _derivation(y) is not S.Zero:
                c, q = p.as_poly(y).primitive()

#.........这里部分代码省略.........
                                terms.add(asin(sqrt(-M[a]/M[b])*x))

                        M = g.base.match(a*x**2 - b)

                        if M is not None and M[b].is_positive:
                            if M[a].is_positive:
                                terms.add(acosh(sqrt(M[a]/M[b])*x))
                            elif M[a].is_negative:
                                terms.add((-M[b]/2*sqrt(-M[a])*\
                                           atan(sqrt(-M[a])*x/sqrt(M[a]*x**2-M[b]))))

        else:
            terms |= set(hints)

    for g in set(terms):
        terms |= components(cancel(g.diff(x)), x)

    V = _symbols('x', len(terms))

    mapping = dict(zip(terms, V))

    rev_mapping = {}

    for k, v in mapping.iteritems():
        rev_mapping[v] = k

    def substitute(expr):
        return expr.subs(mapping)

    diffs = [ substitute(cancel(g.diff(x))) for g in terms ]

    denoms = [ g.as_numer_denom()[1] for g in diffs ]
    try:
        denom = reduce(lambda p, q: lcm(p, q, *V), denoms)
    except PolynomialError:
        # lcm can fail with this. See issue 1418.
        return None

    numers = [ cancel(denom * g) for g in diffs ]

    def derivation(h):
        return Add(*[ d * h.diff(v) for d, v in zip(numers, V) ])

    def deflation(p):
        for y in V:
            if not p.has(y):
                continue

            if derivation(p) is not S.Zero:
                c, q = p.as_poly(y).primitive()
                return deflation(c)*gcd(q, q.diff(y)).as_expr()
        else:
            return p

    def splitter(p):
        for y in V:
            if not p.has(y):
                continue

            if derivation(y) is not S.Zero:
                c, q = p.as_poly(y).primitive()

                q = q.as_expr()

                h = gcd(q, derivation(q), y)
                s = quo(h, gcd(q, q.diff(y), y), y)

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))
       C0*gamma(1 + n) + C1*2**n

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

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

    if f.is_Add:
        F = f.args
    else:
        F = [f]

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

    h_part = {}
    i_part = S.Zero

    for g in F:
        if g.is_Mul:
            G = g.args
        else:
            G = [g]

        coeff = S.One
        kspec = None

        for h in 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:
            if kspec in h_part:
                h_part[kspec] += coeff
            else:
                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):
#.........这里部分代码省略.........

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 parametric_log_deriv_heu(fa, fd, wa, wd, DE, c1=None):
    """
    Parametric logarithmic derivative heuristic.

    Given a derivation D on k[t], f in k(t), and a hyperexponential monomial
    theta over k(t), raises either NotImplementedError, in which case the
    heuristic failed, or returns None, in which case it has proven that no
    solution exists, or returns a solution (n, m, v) of the equation
    n*f == Dv/v + m*Dtheta/theta, with v in k(t)* and n, m in ZZ with n != 0.

    If this heuristic fails, the structure theorem approach will need to be
    used.

    The argument w == Dtheta/theta
    """
    # TODO: finish writing this and write tests
    c1 = c1 or Dummy('c1')

    p, a = fa.div(fd)
    q, b = wa.div(wd)

    B = max(0, derivation(DE.t, DE).degree(DE.t) - 1)
    C = max(p.degree(DE.t), q.degree(DE.t))

    if q.degree(DE.t) > B:
        eqs = [p.nth(i) - c1*q.nth(i) for i in range(B + 1, C + 1)]
        s = solve(eqs, c1)
        if not s or not s[c1].is_Rational:
            # deg(q) > B, no solution for c.
            return None

        N, M = s[c1].as_numer_denom()  # N and M are integers
        N, M = Poly(N, DE.t), Poly(M, DE.t)

        nfmwa = N*fa*wd - M*wa*fd
        nfmwd = fd*wd
        Qv = is_log_deriv_k_t_radical_in_field(N*fa*wd - M*wa*fd, fd*wd, DE,
            'auto')
        if Qv is None:
            # (N*f - M*w) is not the logarithmic derivative of a k(t)-radical.
            return None

        Q, e, v = Qv
        if e != 1:
            return None

        if Q.is_zero or v.is_zero:
            # Q == 0 or v == 0.
            return None

        return (Q*N, Q*M, v)

    if p.degree(DE.t) > B:
        return None

    c = lcm(fd.as_poly(DE.t).LC(), wd.as_poly(DE.t).LC())
    l = fd.monic().lcm(wd.monic())*Poly(c, DE.t)
    ln, ls = splitfactor(l, DE)
    z = ls*ln.gcd(ln.diff(DE.t))

    if not z.has(DE.t):
        raise NotImplementedError("parametric_log_deriv_heu() "
            "heuristic failed: z in k.")

    u1, r1 = (fa*l.quo(fd)).div(z)  # (l*f).div(z)
    u2, r2 = (wa*l.quo(wd)).div(z)  # (l*w).div(z)

    eqs = [r1.nth(i) - c1*r2.nth(i) for i in range(z.degree(DE.t))]
    s = solve(eqs, c1)
    if not s or not s[c1].is_Rational:
        # deg(q) <= B, no solution for c.
        return None

    M, N = s[c1].as_numer_denom()

    nfmwa = N.as_poly(DE.t)*fa*wd - M.as_poly(DE.t)*wa*fd
    nfmwd = fd*wd
    Qv = is_log_deriv_k_t_radical_in_field(nfmwa, nfmwd, DE)
    if Qv is None:
        # (N*f - M*w) is not the logarithmic derivative of a k(t)-radical.
        return None

    Q, v = Qv

    if Q.is_zero or v.is_zero:
        # Q == 0 or v == 0.
        return None

    return (Q*N, Q*M, v)

#.........这里部分代码省略.........
            '''
            newe, r = divmod(common_b[b], b[1])
            if not r:
                common_b.pop(b)
                if newe:
                    for m in Mul.make_args(b[0]**newe):
                        b, e = bkey(m)
                        if b not in common_b:
                            common_b[b] = 0
                        common_b[b] += e
                        if b[1] != 1:
                            bases.append(b)
        # ---------------- end of helper functions

        # assemble a dictionary of the factors having a Rational power
        common_b = {}
        done = []
        bases = []
        for b, e in c_powers:
            b, e = bkey(b, e)
            if b in common_b:
                common_b[b] = common_b[b] + e
            else:
                common_b[b] = e
            if b[1] != 1 and b[0].is_Mul:
                bases.append(b)
        bases.sort(key=default_sort_key)  # this makes tie-breaking canonical
        bases.sort(key=measure, reverse=True)  # handle longest first
        for base in bases:
            if base not in common_b:  # it may have been removed already
                continue
            b, exponent = base
            last = False  # True when no factor of base is a radical
            qlcm = 1  # the lcm of the radical denominators
            while True:
                bstart = b
                qstart = qlcm

                bb = []  # list of factors
                ee = []  # (factor's expo. and it's current value in common_b)
                for bi in Mul.make_args(b):
                    bib, bie = bkey(bi)
                    if bib not in common_b or common_b[bib] < bie:
                        ee = bb = []  # failed
                        break
                    ee.append([bie, common_b[bib]])
                    bb.append(bib)
                if ee:
                    # find the number of integral extractions possible
                    # e.g. [(1, 2), (2, 2)] -> min(2/1, 2/2) -> 1
                    min1 = ee[0][1]//ee[0][0]
                    for i in range(1, len(ee)):
                        rat = ee[i][1]//ee[i][0]
                        if rat < 1:
                            break
                        min1 = min(min1, rat)
                    else:
                        # update base factor counts
                        # e.g. if ee = [(2, 5), (3, 6)] then min1 = 2
                        # and the new base counts will be 5-2*2 and 6-2*3
                        for i in range(len(bb)):
                            common_b[bb[i]] -= min1*ee[i][0]
                            update(bb[i])
                        # update the count of the base
                        # e.g. x**2*y*sqrt(x*sqrt(y)) the count of x*sqrt(y)
                        # will increase by 4 to give bkey (x*sqrt(y), 2, 5)

        def lcm(f, g):
            from sympy.polys import lcm

            return f.__class__(lcm(f.ex, f.__class__(g).ex))

def spoly(f, g):
    """Computes the S-polynomial of f and g"""
    lt_f = _LT(f)
    lt_g = _LT(g)
    common = lcm(lt_f, lt_g)
    return ((common * g) / lt_g - (common * f) / lt_f).expand()