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

    The inhomogeneous part can be either hypergeometric or a sum
    of a fixed number of pairwise dissimilar hypergeometric terms.

    The algorithm performs three basic steps:

        (1) Group together similar hypergeometric terms in the
            inhomogeneous part of `\operatorname{L} y = f`, and find
            particular solution using Abramov's algorithm.

        (2) Compute generating set of `\operatorname{L}` and find basis
            in it, so that all solutions are linearly independent.

        (3) Form final solution with the number of arbitrary
            constants equal to dimension of basis of `\operatorname{L}`.

    Term `a(n)` is hypergeometric if it is annihilated by first order
    linear difference equations with polynomial coefficients or, in
    simpler words, if consecutive term ratio is a rational function.

    The output of this procedure is a linear combination of fixed
    number of hypergeometric terms. However the underlying method
    can generate larger class of solutions - D'Alembertian terms.

    Note also that this method not only computes the kernel of the
    inhomogeneous equation, but also reduces in to a basis so that
    solutions generated by this procedure are linearly independent

    Examples
    ========

    >>> from sympy.solvers import rsolve_hyper
    >>> from sympy.abc import x

    >>> rsolve_hyper([-1, -1, 1], 0, x)
    C0*(1/2 + sqrt(5)/2)**x + C1*(-sqrt(5)/2 + 1/2)**x

    >>> rsolve_hyper([-1, 1], 1 + x, x)
    C0 + x*(x + 1)/2

    References
    ==========

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

    .. [2] M. Petkovsek, H. S. Wilf, D. Zeilberger, A = B, 1996.
Given linear recurrence operator `\operatorname{L}` of order `k`
    with polynomial coefficients and inhomogeneous equation
    `\operatorname{L} y = f` we seek for all hypergeometric solutions
    over field `K` of characteristic zero.
Given linear recurrence operator \ operatorname { L } of order k with polynomial coefficients and inhomogeneous equation \ operatorname { L } y = f we seek for all hypergeometric solutions over field K of characteristic zero .
def rsolve_hyper(coeffs, f, n, **hints):
    """
    Given linear recurrence operator `\operatorname{L}` of order `k`
    with polynomial coefficients and inhomogeneous equation
    `\operatorname{L} y = f` we seek for all hypergeometric solutions
    over field `K` of characteristic zero.

    The inhomogeneous part can be either hypergeometric or a sum
    of a fixed number of pairwise dissimilar hypergeometric terms.

    The algorithm performs three basic steps:

        (1) Group together similar hypergeometric terms in the
            inhomogeneous part of `\operatorname{L} y = f`, and find
            particular solution using Abramov's algorithm.

        (2) Compute generating set of `\operatorname{L}` and find basis
            in it, so that all solutions are linearly independent.

        (3) Form final solution with the number of arbitrary
            constants equal to dimension of basis of `\operatorname{L}`.

    Term `a(n)` is hypergeometric if it is annihilated by first order
    linear difference equations with polynomial coefficients or, in
    simpler words, if consecutive term ratio is a rational function.

    The output of this procedure is a linear combination of fixed
    number of hypergeometric terms. However the underlying method
    can generate larger class of solutions - D'Alembertian terms.

    Note also that this method not only computes the kernel of the
    inhomogeneous equation, but also reduces in to a basis so that
    solutions generated by this procedure are linearly independent

    Examples
    ========

    >>> from sympy.solvers import rsolve_hyper
    >>> from sympy.abc import x

    >>> rsolve_hyper([-1, -1, 1], 0, x)
    C0*(1/2 + sqrt(5)/2)**x + C1*(-sqrt(5)/2 + 1/2)**x

    >>> rsolve_hyper([-1, 1], 1 + x, x)
    C0 + x*(x + 1)/2

    References
    ==========

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

    .. [2] M. Petkovsek, H. S. Wilf, D. Zeilberger, A = B, 1996.
    """
    coeffs = list(map(sympify, coeffs))

    f = sympify(f)

    r, kernel, symbols = len(coeffs) - 1, [], set()

    if not f.is_zero:
        if f.is_Add:
            similar = {}

            for g in f.expand().args:
                if not g.is_hypergeometric(n):
                    return None

                for h in similar.keys():
                    if hypersimilar(g, h, n):
                        similar[h] += g
                        break
                else:
                    similar[g] = S.Zero

            inhomogeneous = []

            for g, h in similar.items():
                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 = Dummy('Z')

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

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

    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

        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).keys():
            if z.is_zero:
                continue

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

            if C is not None and C is not S.Zero:
                symbols |= set(s)

                ratio = z * A * C.subs(n, n + 1) / B / C
                ratio = simplify(ratio)
                # If there is a nonnegative root in the denominator of the ratio,
                # this indicates that the term y(n_root) is zero, and one should
                # start the product with the term y(n_root + 1).
                n0 = 0
                for n_root in roots(ratio.as_numer_denom()[1], n).keys():
                    if (n0 < (n_root + 1)) == True:
                        n0 = n_root + 1
                K = product(ratio, (n, n0, n - 1))
                if K.has(factorial, FallingFactorial, RisingFactorial):
                    K = simplify(K)

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

    kernel.sort(key=default_sort_key)
    sk = list(zip(numbered_symbols('C'), kernel))

    if sk:
        for C, ker in sk:
            result += C * ker
    else:
        return None

    if hints.get('symbols', False):
        symbols |= set([s for s, k in sk])
        return (result, list(symbols))
    else:
        return result
def rsolve_hyper ( coeffs , f , n , * * hints ) : coeffs = list ( map ( sympify , coeffs ) ) f = sympify ( f ) r , kernel , symbols = len ( coeffs ) - 1 , [ ] , set ( ) if not f . is_zero : if f . is_Add : similar = { } for g in f . expand ( ) . args : if not g . is_hypergeometric ( n ) : return None for h in similar . keys ( ) : if hypersimilar ( g , h , n ) : similar [ h ] += g break else : similar [ g ] = S . Zero inhomogeneous = [ ] for g , h in similar . items ( ) : 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 = Dummy ( 'Z' ) p , q = coeffs [ 0 ] , coeffs [ r ] . subs ( n , n - r + 1 ) p_factors = [ z for z in roots ( p , n ) . keys ( ) ] q_factors = [ z for z in roots ( q , n ) . keys ( ) ] 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 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 ) . keys ( ) : if z . is_zero : continue ( C , s ) = rsolve_poly ( [ polys [ i ] * z ** i for i in xrange ( r + 1 ) ] , 0 , n , symbols = True ) if C is not None and C is not S . Zero : symbols |= set ( s ) ratio = z * A * C . subs ( n , n + 1 ) / B / C ratio = simplify ( ratio ) # If there is a nonnegative root in the denominator of the ratio, # this indicates that the term y(n_root) is zero, and one should # start the product with the term y(n_root + 1). n0 = 0 for n_root in roots ( ratio . as_numer_denom ( ) [ 1 ] , n ) . keys ( ) : if ( n0 < ( n_root + 1 ) ) == True : n0 = n_root + 1 K = product ( ratio , ( n , n0 , n - 1 ) ) if K . has ( factorial , FallingFactorial , RisingFactorial ) : K = simplify ( K ) if casoratian ( kernel + [ K ] , n , zero = False ) != 0 : kernel . append ( K ) kernel . sort ( key = default_sort_key ) sk = list ( zip ( numbered_symbols ( 'C' ) , kernel ) ) if sk : for C , ker in sk : result += C * ker else : return None if hints . get ( 'symbols' , False ) : symbols |= set ( [ s for s , k in sk ] ) return ( result , list ( symbols ) ) else : return result