def decompose_clifford_tableau_to_operations(
    qubits: list[cirq.Qid], clifford_tableau: qis.CliffordTableau
) -> list[ops.Operation]:
    """Decompose an n-qubit Clifford Tableau into a list of one/two qubit operations.

    The implementation is based on Theorem 8 in [1].
    [1] S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
        Phys. Rev. A 70, 052328 (2004). https://arxiv.org/abs/quant-ph/0406196

    Args:
        qubits: The list of qubits being operated on.
        clifford_tableau: The Clifford Tableau for decomposition.

    Returns:
        A list of operations reconstructs the same Clifford tableau.

    Raises:
        ValueError: The length of input qubit mismatch with the size of tableau.
    """
    if len(qubits) != clifford_tableau.n:
        raise ValueError("The number of qubits must be the same as the number of Clifford Tableau.")
    assert (
        clifford_tableau._validate()
    ), "The provided clifford_tableau must satisfy the symplectic property."

    t: qis.CliffordTableau = clifford_tableau.copy()
    operations: list[ops.Operation] = []
    args = sim.CliffordTableauSimulationState(
        tableau=t, qubits=qubits, prng=np.random.RandomState()
    )

    _X_with_ops = functools.partial(_X, args=args, operations=operations, qubits=qubits)
    _Z_with_ops = functools.partial(_Z, args=args, operations=operations, qubits=qubits)
    _H_with_ops = functools.partial(_H, args=args, operations=operations, qubits=qubits)
    _S_with_ops = functools.partial(_Sdg, args=args, operations=operations, qubits=qubits)
    _CNOT_with_ops = functools.partial(_CNOT, args=args, operations=operations, qubits=qubits)
    _SWAP_with_ops = functools.partial(_SWAP, args=args, operations=operations, qubits=qubits)

    # The procedure is based on Theorem 8 in
    # [1] S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
    #     Phys. Rev. A 70, 052328 (2004). https://arxiv.org/abs/quant-ph/0406196
    # with modification by doing it row-by-row instead.

    # Suppose we have a Clifford Tableau:
    #                   Xs  Zs
    # Destabilizers:  [ A | B ]
    # Stabilizers:    [ C | D ]
    for i in range(t.n):
        # Step 1a: Make the diagonal element of A equal to 1 by Hadamard gate if necessary.
        if not t.xs[i, i] and t.zs[i, i]:
            _H_with_ops(i)
        # Step 1b: Make the diagonal element of A equal to 1 by SWAP gate if necessary.
        if not t.xs[i, i]:
            for j in range(i + 1, t.n):
                if t.xs[i, j]:
                    _SWAP_with_ops(i, j)
                    break
        # Step 1c: We may still not be able to find non-zero element in whole Xs row. Then,
        # apply swap + Hadamard from zs. It is guaranteed to find one by lemma 5 in [1].
        if not t.xs[i, i]:
            for j in range(i + 1, t.n):
                if t.zs[i, j]:
                    _H_with_ops(j)
                    _SWAP_with_ops(i, j)
                    break

        # Step 2: Eliminate the elements in A By CNOT and phase gate (i-th row)
        # first i rows of destabilizers: [ I  0 | 0  0 ]
        _ = [_CNOT_with_ops(i, j) for j in range(i + 1, t.n) if t.xs[i, j]]
        if np.any(t.zs[i, i:]):
            if not t.zs[i, i]:
                _S_with_ops(i)
            _ = [_CNOT_with_ops(j, i) for j in range(i + 1, t.n) if t.zs[i, j]]
            _S_with_ops(i)

        # Step 3: Eliminate the elements in D By CNOT and phase gate (i-th row)
        # first i rows of stabilizers: [ 0  0 | I  0 ]
        _ = [_CNOT_with_ops(j, i) for j in range(i + 1, t.n) if t.zs[i + t.n, j]]
        if np.any(t.xs[i + t.n, i:]):
            # Swap xs and zs
            _H_with_ops(i)
            _ = [_CNOT_with_ops(i, j) for j in range(i + 1, t.n) if t.xs[i + t.n, j]]
            if t.zs[i + t.n, i]:
                _S_with_ops(i)
            _H_with_ops(i)

    # Step 4: Correct the phase of tableau
    _ = [_Z_with_ops(i) for i, p in enumerate(t.rs[: t.n]) if p]
    _ = [_X_with_ops(i) for i, p in enumerate(t.rs[t.n :]) if p]

    # Step 5: invert the operations by reversing the order: (AB)^{+} = B^{+} A^{+}.
    # Note only S gate is not self-adjoint.
    return operations[::-1]
Decompose an n-qubit Clifford Tableau into a list of one/two qubit operations.

    The implementation is based on Theorem 8 in [1].
    [1] S. Aaronson, D. Gottesman, *Improved Simulation of Stabilizer Circuits*,
        Phys. Rev. A 70, 052328 (2004). https://arxiv.org/abs/quant-ph/0406196

    Args:
        qubits: The list of qubits being operated on.
        clifford_tableau: The Clifford Tableau for decomposition.

    Returns:
        A list of operations reconstructs the same Clifford tableau.

    Raises:
        ValueError: The length of input qubit mismatch with the size of tableau.
    
cirq-core/cirq/transformers/analytical_decompositions/clifford_decomposition.py