code: def poly_linear_constraints(p, d):
      """
      Given p = [p1, ..., pm] in k[t]^m and d in k[t], return
      q = [q1, ..., qm] in k[t]^m and a matrix M with entries in k such
      that Sum(ci*pi, (i, 1, m)), for c1, ..., cm in k, is divisible
      by d if and only if (c1, ..., cm) is a solution of Mx = 0, in
      which case the quotient is Sum(ci*qi, (i, 1, m)).
      """
      m = len(p)
      q, r = zip(*[pi.div(d) for pi in p])
  
      if not all(ri.is_zero for ri in r):
          n = max(ri.degree() for ri in r)
          M = Matrix(n + 1, m, lambda i, j: r[j].nth(i), d.gens)
      else:
          M = Matrix(0, m, [], d.gens)  # No constraints.
  
      return q, M
  
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 5 for calling a global function with example usage: poly_linear_constraints(p, d) and returns: q,  M