code: def param_poly_rischDE(a, b, q, n, DE):
      """Polynomial solutions of a parametric Risch differential equation.
  
      Explanation
      ===========
  
      Given a derivation D in k[t], a, b in k[t] relatively prime, and q
      = [q1, ..., qm] in k[t]^m, return h = [h1, ..., hr] in k[t]^r and
      a matrix A with m + r columns and entries in Const(k) such that
      a*Dp + b*p = Sum(ci*qi, (i, 1, m)) has a solution p of degree <= n
      in k[t] with c1, ..., cm in Const(k) if and only if p = Sum(dj*hj,
      (j, 1, r)) where d1, ..., dr are in Const(k) and (c1, ..., cm,
      d1, ..., dr) is a solution of Ax == 0.
      """
      m = len(q)
      if n < 0:
          # Only the trivial zero solution is possible.
          # Find relations between the qi.
          if all(qi.is_zero for qi in q):
              return [], zeros(1, m, DE.t)  # No constraints.
  
          N = max([qi.degree(DE.t) for qi in q])
          M = Matrix(N + 1, m, lambda i, j: q[j].nth(i), DE.t)
          A, _ = constant_system(M, zeros(M.rows, 1, DE.t), DE)
  
          return [], A
  
      if a.is_ground:
          # Normalization: a = 1.
          a = a.LC()
          b, q = b.quo_ground(a), [qi.quo_ground(a) for qi in q]
  
          if not b.is_zero and (
              DE.case == "base"
              or b.degree() > max(0, DE.d.degree() - 1)
          ):
              return prde_no_cancel_b_large(b, q, n, DE)
  
          elif (b.is_zero or b.degree() < DE.d.degree() - 1) and (
              DE.case == "base" or DE.d.degree() >= 2
          ):
              return prde_no_cancel_b_small(b, q, n, DE)
  
          elif (
              DE.d.degree() >= 2
              and b.degree() == DE.d.degree() - 1
              and n > -b.as_poly().LC() / DE.d.as_poly().LC()
          ):
              raise NotImplementedError(
                  "prde_no_cancel_b_equal() is " "not yet implemented."
              )
  
          else:
              # Liouvillian cases
              if DE.case in ("primitive", "exp"):
                  return prde_cancel_liouvillian(b, q, n, DE)
              else:
                  raise NotImplementedError(
                      "non-linear and hypertangent "
                      "cases have not yet been implemented"
                  )
  
      # else: deg(a) > 0
  
      # Iterate SPDE as long as possible cumulating coefficient
      # and terms for the recovery of original solutions.
      alpha, beta = a.one, [a.zero] * m
      while n >= 0:  # and a, b relatively prime
          a, b, q, r, n = prde_spde(a, b, q, n, DE)
          beta = [betai + alpha * ri for betai, ri in zip(beta, r)]
          alpha *= a
          # Solutions p of a*Dp + b*p = Sum(ci*qi) correspond to
          # solutions alpha*p + Sum(ci*betai) of the initial equation.
          d = a.gcd(b)
          if not d.is_ground:
              break
  
      # a*Dp + b*p = Sum(ci*qi) may have a polynomial solution
      # only if the sum is divisible by d.
  
      qq, M = poly_linear_constraints(q, d)
      # qq = [qq1, ..., qqm] where qqi = qi.quo(d).
      # M is a matrix with m columns an entries in k.
      # Sum(fi*qi, (i, 1, m)), where f1, ..., fm are elements of k, is
      # divisible by d if and only if M*Matrix([f1, ..., fm]) == 0,
      # in which case the quotient is Sum(fi*qqi).
  
      A, _ = constant_system(M, zeros(M.rows, 1, DE.t), DE)
      # A is a matrix with m columns and entries in Const(k).
      # Sum(ci*qqi) is Sum(ci*qi).quo(d), and the remainder is zero
      # for c1, ..., cm in Const(k) if and only if
      # A*Matrix([c1, ...,cm]) == 0.
  
      V = A.nullspace()
      # V = [v1, ..., vu] where each vj is a column matrix with
      # entries aj1, ..., ajm in Const(k).
      # Sum(aji*qi) is divisible by d with exact quotient Sum(aji*qqi).
      # Sum(ci*qi) is divisible by d if and only if ci = Sum(dj*aji)
      # (i = 1, ..., m) for some d1, ..., du in Const(k).
      # In that case, solutions of
      #     a*Dp + b*p = Sum(ci*qi) = Sum(dj*Sum(aji*qi))
      # are the same as those of
      #     (a/d)*Dp + (b/d)*p = Sum(dj*rj)
      # where rj = Sum(aji*qqi).
  
      if not V:  # No non-trivial solution.
          return [], eye(m, DE.t)  # Could return A, but this has
          # the minimum number of rows.
  
      Mqq = Matrix([qq])  # A single row.
      r = [(Mqq * vj)[0] for vj in V]  # [r1, ..., ru]
  
      # Solutions of (a/d)*Dp + (b/d)*p = Sum(dj*rj) correspond to
      # solutions alpha*p + Sum(Sum(dj*aji)*betai) of the initial
      # equation. These are equal to alpha*p + Sum(dj*fj) where
      # fj = Sum(aji*betai).
      Mbeta = Matrix([beta])
      f = [(Mbeta * vj)[0] for vj in V]  # [f1, ..., fu]
  
      #
      # Solve the reduced equation recursively.
      #
      g, B = param_poly_rischDE(a.quo(d), b.quo(d), r, n, DE)
  
      # g = [g1, ..., gv] in k[t]^v and and B is a matrix with u + v
      # columns and entries in Const(k) such that
      # (a/d)*Dp + (b/d)*p = Sum(dj*rj) has a solution p of degree <= n
      # in k[t] if and only if p = Sum(ek*gk) where e1, ..., ev are in
      # Const(k) and B*Matrix([d1, ..., du, e1, ..., ev]) == 0.
      # The solutions of the original equation are then
      # Sum(dj*fj, (j, 1, u)) + alpha*Sum(ek*gk, (k, 1, v)).
  
      # Collect solution components.
      h = f + [alpha * gk for gk in g]
  
      # Build combined relation matrix.
      A = -eye(m, DE.t)
      for vj in V:
          A = A.row_join(vj)
      A = A.row_join(zeros(m, len(g), DE.t))
      A = A.col_join(zeros(B.rows, m, DE.t).row_join(B))
  
      return h, A
  
imports: import itertools
  from functools import reduce
  from sympy.core.intfunc import ilcm
  from sympy.core import Dummy, Add, Mul, Pow, S
  from sympy.integrals.rde import order_at, order_at_oo, weak_normalizer, bound_degree
  from sympy.integrals.risch import gcdex_diophantine, frac_in, derivation, residue_reduce, splitfactor, residue_reduce_derivation, DecrementLevel, recognize_log_derivative
  from sympy.polys import Poly, lcm, cancel, sqf_list
  from sympy.polys.polymatrix import PolyMatrix
  from sympy.solvers import solve
label_desc: Use image node_id 11 for calling a global function with example usage: param_poly_rischDE(a, b, q, n, DE) and returns: h,  A, list,  A, list,  eye, list,  zeros, prde_no_cancel_b_large, prde_no_cancel_b_small, prde_cancel_liouvillian