code: def _hermite_normal_form_modulo_D(A, D):
      r"""
      Perform the mod *D* Hermite Normal Form reduction algorithm on
      :py:class:`~.DomainMatrix` *A*.
  
      Explanation
      ===========
  
      If *A* is an $m \times n$ matrix of rank $m$, having Hermite Normal Form
      $W$, and if *D* is any positive integer known in advance to be a multiple
      of $\det(W)$, then the HNF of *A* can be computed by an algorithm that
      works mod *D* in order to prevent coefficient explosion.
  
      Parameters
      ==========
  
      A : :py:class:`~.DomainMatrix` over :ref:`ZZ`
          $m \times n$ matrix, having rank $m$.
      D : :ref:`ZZ`
          Positive integer, known to be a multiple of the determinant of the
          HNF of *A*.
  
      Returns
      =======
  
      :py:class:`~.DomainMatrix`
          The HNF of matrix *A*.
  
      Raises
      ======
  
      DMDomainError
          If the domain of the matrix is not :ref:`ZZ`, or
          if *D* is given but is not in :ref:`ZZ`.
  
      DMShapeError
          If the matrix has more rows than columns.
  
      References
      ==========
  
      .. [1] Cohen, H. *A Course in Computational Algebraic Number Theory.*
         (See Algorithm 2.4.8.)
  
      """
      if not A.domain.is_ZZ:
          raise DMDomainError("Matrix must be over domain ZZ.")
      if not ZZ.of_type(D) or D < 1:
          raise DMDomainError(
              "Modulus D must be positive element of domain ZZ."
          )
  
      def add_columns_mod_R(m, R, i, j, a, b, c, d):
          # replace m[:, i] by (a*m[:, i] + b*m[:, j]) % R
          # and m[:, j] by (c*m[:, i] + d*m[:, j]) % R
          for k in range(len(m)):
              e = m[k][i]
              m[k][i] = symmetric_residue((a * e + b * m[k][j]) % R, R)
              m[k][j] = symmetric_residue((c * e + d * m[k][j]) % R, R)
  
      W = defaultdict(dict)
  
      m, n = A.shape
      if n < m:
          raise DMShapeError(
              "Matrix must have at least as many columns as rows."
          )
      A = A.to_dense().rep.to_ddm().copy()
      k = n
      R = D
      for i in range(m - 1, -1, -1):
          k -= 1
          for j in range(k - 1, -1, -1):
              if A[i][j] != 0:
                  u, v, d = _gcdex(A[i][k], A[i][j])
                  r, s = A[i][k] // d, A[i][j] // d
                  add_columns_mod_R(A, R, k, j, u, v, -s, r)
          b = A[i][k]
          if b == 0:
              A[i][k] = b = R
          u, v, d = _gcdex(b, R)
          for ii in range(m):
              W[ii][i] = u * A[ii][k] % R
          if W[i][i] == 0:
              W[i][i] = R
          for j in range(i + 1, m):
              q = W[i][j] // W[i][i]
              add_columns(W, j, i, 1, -q, 0, 1)
          R //= d
      return DomainMatrix(W, (m, m), ZZ).to_dense()
  
file_path: sympy/sympy/polys/matrices/normalforms.py
imports: from collections import defaultdict
  from .domainmatrix import DomainMatrix
  from .exceptions import DMDomainError, DMShapeError
  from sympy.ntheory.modular import symmetric_residue
  from sympy.polys.domains import QQ, ZZ
label_desc: Use image node_id 6 for calling a global function with example usage: _hermite_normal_form_modulo_D(A, D) and returns: DomainMatrix