code: def rsolve_hyper(coeffs, f, n, **hints):
      r"""
      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*(1/2 + sqrt(5)/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 = [g + h for g, h in similar.items()]
          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 range(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 range(r + 1):
                  polys[j] *= Mul(*(denoms[:j] + denoms[j + 1 :]))
  
              # FIXME: The call to rsolve_ratio below should suffice (rsolve_poly
              # call can be removed) but the XFAIL test_rsolve_ratio_missed must
              # be fixed first.
              R = rsolve_ratio(polys, Mul(*denoms), n, symbols=True)
              if R is not None:
                  R, syms = R
                  if syms:
                      R = R.subs(zip(syms, [0] * len(syms)))
              else:
                  R = rsolve_poly(polys, Mul(*denoms), n)
  
              if R:
                  inhomogeneous[i] *= R
              else:
                  return None
  
              result = Add(*inhomogeneous)
              result = simplify(result)
      else:
          result = S.Zero
  
      Z = Dummy("Z")
  
      p, q = coeffs[0], coeffs[r].subs(n, n - r + 1)
  
      p_factors = list(roots(p, n).keys())
      q_factors = list(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 range(r + 1):
              a = Mul(*[A.subs(n, n + j) for j in range(i)])
              b = Mul(*[B.subs(n, n + j) for j in range(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())
  
          if degrees:
              d, poly = max(degrees), S.Zero
          else:
              return None
  
          for i in range(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
  
              recurr_coeffs = [
                  polys[i].as_expr() * z**i for i in range(r + 1)
              ]
              if d == 0 and 0 != Add(
                  *[recurr_coeffs[j] * j for j in range(1, r + 1)]
              ):
                  # faster inline check (than calling rsolve_poly) for a
                  # constant solution to a constant coefficient recurrence.
                  sol = [Symbol("C" + str(len(symbols)))]
              else:
                  sol, syms = rsolve_poly(
                      recurr_coeffs, 0, n, len(symbols), symbols=True
                  )
                  sol = sol.collect(syms)
                  sol = [sol.coeff(s) for s in syms]
  
              for C in sol:
                  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 n_root.has(I):
                          return None
                      elif (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))
  
      for C, ker in sk:
          result += C * ker
  
      if hints.get("symbols", False):
          # XXX: This returns the symbols in a non-deterministic order
          symbols |= {s for s, k in sk}
          return (result, list(symbols))
      else:
          return result
  
imports: from collections import defaultdict
  from sympy.concrete import product
  from sympy.core.singleton import S
  from sympy.core.numbers import Rational, I
  from sympy.core.symbol import Symbol, Wild, Dummy
  from sympy.core.relational import Equality
  from sympy.core.add import Add
  from sympy.core.mul import Mul
  from sympy.core.sorting import default_sort_key
  from sympy.core.sympify import sympify
  from sympy.simplify import simplify, hypersimp, hypersimilar
  from sympy.solvers import solve, solve_undetermined_coeffs
  from sympy.polys import Poly, quo, gcd, lcm, roots, resultant
  from sympy.functions import binomial, factorial, FallingFactorial, RisingFactorial
  from sympy.matrices import Matrix, casoratian
  from sympy.utilities.iterables import numbered_symbols
label_desc: Use image node_id 3 for calling a global function with example usage: rsolve_hyper(coeffs, f, n) and returns: result,  list, result, None, None, None, None, None