def _decomp_3sqrt_iswap_matrices(
    kak: cirq.KakDecomposition, atol: float = 1e-8
) -> tuple[Sequence[tuple[np.ndarray, np.ndarray]], complex]:
    """Returns the single-qubit matrices for the 3-SQRT_ISWAP decomposition.

    Assumes any canonical x, y, z.  Three sqrt-iSWAP gates are only needed if
    x < y + |z|.  Only two are needed for other gates (most cases).

    References:
        Towards ultra-high fidelity quantum operations: SQiSW gate as a native
        two-qubit gate
        https://arxiv.org/abs/2105.06074
    """
    # This somewhat follows the else-branch of procedure DECOMP(U) in Algorithm 1 of the paper.
    # However the canonicalization conditions are different from the paper and allow any Weyl
    # coordinate to be synthesized with 3 sqrt-iSWAPs.
    #
    # This method breaks the 3-sqrt-iSWAP synthesis problem into a sum of the 1-sqrt-iSWAP and
    # 2-sqrt-iSWAP problems.  This works because given two 2-qubit unitaries, U and V, with (not
    # necessarily canonical) Weyl coordinates (x1, y1, z1) and (x2, y2, z2), both products U*V and
    # V*U will have (non-canonical) Weyl coordinates (x1+x2, y1+y2, z1+z2) if both U and V are
    # diagonal in the magic basis (i.e. all single-qubit operations of the pre-canonicalized KAK
    # decomposition are identity).
    x, y, z = kak.interaction_coefficients
    b0, b1 = kak.single_qubit_operations_before
    a0, a1 = kak.single_qubit_operations_after

    # Find x1, y1, z1, x2, y2, z2
    # such that x1+x2=x, y1+y2=y, z1+z2=z
    # where x1, y1, z1 are implementable by one sqrt-iSWAP gate
    # and x2, y2, z2 implementable by two sqrt-iSWAP gates
    # No error tolerance needed
    ieq1 = y > np.pi / 8
    ieq2 = z < 0
    if ieq1:
        if ieq2:
            # Non-canonical Weyl coordinates for the single sqrt-iSWAP
            x1, y1, z1 = 0.0, np.pi / 8, -np.pi / 8
        else:
            x1, y1, z1 = 0.0, np.pi / 8, np.pi / 8
    else:
        x1, y1, z1 = -np.pi / 8, np.pi / 8, 0.0
    # Non-canonical Weyl coordinates for the two sqrt-iSWAP decomposition
    x2, y2, z2 = x - x1, y - y1, z - z1

    # Find fixup single-qubit gates for the canonical (i.e. diagonal in the magic basis)
    # decompositions
    kak1 = linalg.kak_canonicalize_vector(x1, y1, z1, atol)
    kak2 = linalg.kak_canonicalize_vector(x2, y2, z2, atol)

    # Compute sub-decompositions
    # F0 and F1 from Algorithm 1 of the paper are not needed
    ((h0, h1), (g0, g1)), phase1 = _decomp_1sqrt_iswap_matrices(kak1, atol)
    ((e0, e1), (c0, c1), (d0, d1)), phase2 = _decomp_2sqrt_iswap_matrices(kak2, atol)

    # There are two valid solutions at this point: kak1 before kak2 or kak2 before kak1
    # Arbitrarily pick kak1 before kak2
    return [  # Pairs of single-qubit unitaries, SQRT_ISWAP between each is implied
        (h0 @ b0, h1 @ b1),
        (e0 @ g0, e1 @ g1),
        (c0, c1),
        (a0 @ d0, a1 @ d1),
    ], kak.global_phase * phase1 * phase2
Returns the single-qubit matrices for the 3-SQRT_ISWAP decomposition.

    Assumes any canonical x, y, z.  Three sqrt-iSWAP gates are only needed if
    x < y + |z|.  Only two are needed for other gates (most cases).

    References:
        Towards ultra-high fidelity quantum operations: SQiSW gate as a native
        two-qubit gate
        https://arxiv.org/abs/2105.06074
    
cirq-core/cirq/transformers/analytical_decompositions/two_qubit_to_sqrt_iswap.py