docs: def diagonalize_real_symmetric_and_sorted_diagonal_matrices(
          symmetric_matrix: np.ndarray,
          diagonal_matrix: np.ndarray,
          *,
          rtol: float = 1e-5,
          atol: float = 1e-8,
          check_preconditions: bool = True) -> np.ndarray:
      """Returns an orthogonal matrix that diagonalizes both given matrices.
  
      The given matrices must commute.
      Guarantees that the sorted diagonal matrix is not permuted by the
      diagonalization (except for nearly-equal values).
  
      Args:
          symmetric_matrix: A real symmetric matrix.
          diagonal_matrix: A real diagonal matrix with entries along the diagonal
              sorted into descending order.
          rtol: Relative numeric error threshold.
          atol: Absolute numeric error threshold.
          check_preconditions: If set, verifies that the input matrices commute
              and are respectively symmetric and diagonal descending.
  
      Returns:
          An orthogonal matrix P such that P.T @ symmetric_matrix @ P is diagonal
          and P.T @ diagonal_matrix @ P = diagonal_matrix (up to tolerance).
  
      Raises:
          ValueError: Matrices don't meet preconditions (e.g. not symmetric).
      """
  
      # Verify preconditions.
      if check_preconditions:
          if (np.any(np.imag(symmetric_matrix)) or
                  not predicates.is_hermitian(symmetric_matrix,
                                              rtol=rtol,
                                              atol=atol)):
              raise ValueError('symmetric_matrix must be real symmetric.')
          if (not predicates.is_diagonal(diagonal_matrix, atol=atol) or
                  np.any(np.imag(diagonal_matrix)) or
                  np.any(diagonal_matrix[:-1, :-1] < diagonal_matrix[1:, 1:])):
              raise ValueError(
                  'diagonal_matrix must be real diagonal descending.')
          if not predicates.commutes(diagonal_matrix,
                                     symmetric_matrix,
                                     rtol=rtol,
                                     atol=atol):
              raise ValueError('Given matrices must commute.')
  
      def similar_singular(i, j):
          return np.allclose(diagonal_matrix[i, i],
                              diagonal_matrix[j, j], rtol=rtol)
  
      # Because the symmetric matrix commutes with the diagonal singulars matrix,
      # the symmetric matrix should be block-diagonal with a block boundary
      # wherever the singular values happen change. So we can use the singular
      # values to extract blocks that can be independently diagonalized.
      ranges = _contiguous_groups(diagonal_matrix.shape[0], similar_singular)
  
      # Build the overall diagonalization by diagonalizing each block.
      p = np.zeros(symmetric_matrix.shape, dtype=np.float64)
      for start, end in ranges:
          block = symmetric_matrix[start:end, start:end]
          p[start:end, start:end] = diagonalize_real_symmetric_matrix(
              block, rtol=rtol, atol=atol)
  
      return p