def _multiplexed_cossin(
    cossin_qubits: list[cirq.Qid], angles: list[float], rot_func: Callable = ops.ry
) -> Iterable[cirq.Operation]:
    """Performs a multiplexed rotation over all qubits in this unitary matrix,

    Uses ry and rz multiplexing for quantum shannon decomposition

    Args:
        cossin_qubits: Subset of total qubits involved in this unitary gate
        angles: List of angles to be multiplexed over for the given type of rotation
        rot_func: Rotation function used for this multiplexing implementation
                    (cirq.ry or cirq.rz)

    Calls:
        No major calls

    Yields: Single operation from OP TREE from set 1- and 2-qubit gates: {ry,rz,CNOT}
    """
    # Most significant qubit is main qubit with rotation function applied
    main_qubit = cossin_qubits[0]

    # All other qubits are control qubits
    control_qubits = cossin_qubits[1:]

    for j in range(len(angles)):
        # The rotation includes a factor of (-1) for each bit in the Gray Code
        #   if the position of that bit is also 1
        # The number of factors of -1 is counted using the 1s in the
        #   binary representation of the (gray(j) & i)
        # Here, i gives the index for the angle, and
        #   j is the iteration of the decomposition
        rotation = sum(
            -angle if bin(_nth_gray(j) & i).count('1') % 2 else angle
            for i, angle in enumerate(angles)
        )

        # Divide by a factor of 2 for each additional select qubit
        # This is due to the halving in the decomposition applied recursively
        rotation = rotation * 2 / len(angles)

        # The XOR of the this gray code with the next will give the 1 for the bit
        #   corresponding to the CNOT select, else 0
        select_string = _nth_gray(j) ^ _nth_gray(j + 1)

        # Find the index number where the bit is 1
        select_qubit = next(i for i in range(len(angles)) if (select_string >> i & 1))

        # Negate the value, since we must index starting at most significant qubit
        # Also the final value will overflow, and it should be the MSB,
        #   so introduce max function
        select_qubit = max(-select_qubit - 1, -len(control_qubits))

        if np.abs(rotation) > 1e-9:
            # Add a rotation on the main qubit
            yield rot_func(rotation).on(main_qubit)

        # Add a CNOT from the select qubit to the main qubit
        yield ops.CNOT(control_qubits[select_qubit], main_qubit)
Performs a multiplexed rotation over all qubits in this unitary matrix,

    Uses ry and rz multiplexing for quantum shannon decomposition

    Args:
        cossin_qubits: Subset of total qubits involved in this unitary gate
        angles: List of angles to be multiplexed over for the given type of rotation
        rot_func: Rotation function used for this multiplexing implementation
                    (cirq.ry or cirq.rz)

    Calls:
        No major calls

    Yields: Single operation from OP TREE from set 1- and 2-qubit gates: {ry,rz,CNOT}
    
cirq-core/cirq/transformers/analytical_decompositions/quantum_shannon_decomposition.py