code: def _diagonalize(M, reals_only=False, sort=False, normalize=False):
      """
      Return (P, D), where D is diagonal and
  
          D = P^-1 * M * P
  
      where M is current matrix.
  
      Parameters
      ==========
  
      reals_only : bool. Whether to throw an error if complex numbers are need
                      to diagonalize. (Default: False)
  
      sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
  
      normalize : bool. If True, normalize the columns of P. (Default: False)
  
      Examples
      ========
  
      >>> from sympy import Matrix
      >>> M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
      >>> M
      Matrix([
      [1,  2, 0],
      [0,  3, 0],
      [2, -4, 2]])
      >>> (P, D) = M.diagonalize()
      >>> D
      Matrix([
      [1, 0, 0],
      [0, 2, 0],
      [0, 0, 3]])
      >>> P
      Matrix([
      [-1, 0, -1],
      [ 0, 0, -1],
      [ 2, 1,  2]])
      >>> P.inv() * M * P
      Matrix([
      [1, 0, 0],
      [0, 2, 0],
      [0, 0, 3]])
  
      See Also
      ========
  
      sympy.matrices.matrixbase.MatrixBase.is_diagonal
      is_diagonalizable
      """
  
      if not M.is_square:
          raise NonSquareMatrixError()
  
      is_diagonalizable, eigenvecs = _is_diagonalizable_with_eigen(
          M, reals_only=reals_only
      )
  
      if not is_diagonalizable:
          raise MatrixError("Matrix is not diagonalizable")
  
      if sort:
          eigenvecs = sorted(eigenvecs, key=default_sort_key)
  
      p_cols, diag = [], []
  
      for val, mult, basis in eigenvecs:
          diag += [val] * mult
          p_cols += basis
  
      if normalize:
          p_cols = [v / v.norm() for v in p_cols]
  
      return M.hstack(*p_cols), M.diag(*diag)
  
imports: from types import FunctionType
  from collections import Counter
  from mpmath import mp, workprec
  from mpmath.libmp.libmpf import prec_to_dps
  from sympy.core.sorting import default_sort_key
  from sympy.core.evalf import DEFAULT_MAXPREC, PrecisionExhausted
  from sympy.core.logic import fuzzy_and, fuzzy_or
  from sympy.core.numbers import Float
  from sympy.core.sympify import _sympify
  from sympy.functions.elementary.miscellaneous import sqrt
  from sympy.polys import roots, CRootOf, ZZ, QQ, EX
  from sympy.polys.matrices import DomainMatrix
  from sympy.polys.matrices.eigen import dom_eigenvects, dom_eigenvects_to_sympy
  from sympy.polys.polytools import gcd
  from .exceptions import MatrixError, NonSquareMatrixError
  from .determinant import _find_reasonable_pivot
  from .utilities import _iszero, _simplify
label_desc: Use image node_id 18 for calling a global function with example usage: _diagonalize(M, reals_only, sort, normalize) and returns: M,  M