def bidiagonalize_unitary_with_special_orthogonals(
        mat: np.ndarray,
        *,
        rtol: float = 1e-5,
        atol: float = 1e-8,
        check_preconditions: bool = True
) -> Tuple[np.ndarray, np.array, np.ndarray]:
    """Finds orthogonal matrices L, R such that L @ matrix @ R is diagonal.

    Args:
        mat: A unitary matrix.
        rtol: Relative numeric error threshold.
        atol: Absolute numeric error threshold.
        check_preconditions: If set, verifies that the input is a unitary matrix
            (to the given tolerances). Defaults to set.

    Returns:
        A triplet (L, d, R) such that L @ mat @ R = diag(d). Both L and R will
        be orthogonal matrices with determinant equal to 1.

    Raises:
        ValueError: Matrices don't meet preconditions (e.g. not real).
    """

    if check_preconditions:
        if not predicates.is_unitary(mat, rtol=rtol, atol=atol):
            raise ValueError('matrix must be unitary.')

    # Note: Because mat is unitary, setting A = real(mat) and B = imag(mat)
    # guarantees that both A @ B.T and A.T @ B are Hermitian.
    left, right = bidiagonalize_real_matrix_pair_with_symmetric_products(
        np.real(mat),
        np.imag(mat),
        rtol=rtol,
        atol=atol,
        check_preconditions=check_preconditions)

    # Convert to special orthogonal w/o breaking diagonalization.
    if np.linalg.det(left) < 0:
        left[0, :] *= -1
    if np.linalg.det(right) < 0:
        right[:, 0] *= -1

    diag = combinators.dot(left, mat, right)

    return left, np.diag(diag), right