code: def _eig_complex_symmetric(M: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
      
      if not np.allclose(M, M.transpose()):
          raise np.linalg.LinAlgError('Not a symmetric matrix')
  
      # The matrix of eigenvectors should be orthogonal.
      # But the standard 'eig' method will fail to return an orthogonal
      # eigenvector matrix when the eigenvalues are degenerate. However,
      # both the real and
      # imaginary part of M must be symmetric with the same orthogonal
      # matrix of eigenvectors. But either the real or imaginary part could
      # vanish. So we use a randomized algorithm where we diagonalize a
      # random linear combination of real and imaginary parts to find the
      # eigenvectors, taking advantage of the 'eigh' subroutine for
      # diagonalizing symmetric matrices.
      # This can fail if we're very unlucky with our random coefficient, so we
      # give the algorithm a few chances to succeed.
  
      # Empirically, never seems to fail on randomly sampled complex
      # symmetric 4x4 matrices.
      # If failure rate is less than 1 in a million, then 16 rounds
      # will have overall failure rate less than 1 in a googol.
      # However, cannot (yet) guarantee that there aren't special cases
      # which have much higher failure rates.
  
      # GEC 2018
  
      max_attempts = 16
      for _ in range(max_attempts):
          c = np.random.uniform(0, 1)
          matrix = c * M.real + (1-c) * M.imag
          _, eigvecs = np.linalg.eigh(matrix)
          eigvecs = np.array(eigvecs, dtype=complex)
          eigvals = np.diag(eigvecs.transpose() @ M @ eigvecs)
  
          # Finish if we got a correct answer.
          reconstructed = eigvecs @ np.diag(eigvals) @ eigvecs.transpose()
          if np.allclose(M, reconstructed):
              return eigvals, eigvecs
  
      # Should never happen. Hopefully.
      raise np.linalg.LinAlgError(
          'Cannot diagonalize complex symmetric matrix.')
docs: Diagonalize a complex symmetric  matrix. The eigenvalues are
      complex, and the eigenvectors form an orthogonal matrix.
  
      Returns:
          eigenvalues, eigenvectors