def then(self, second: CliffordTableau) -> CliffordTableau:
        """Returns a composed CliffordTableau of this tableau and the second tableau.

        Then composed tableau is equal to (up to global phase) the composed
        unitary operation of the two tableaux, i.e. equivalent to applying the unitary
        operation of this CliffordTableau then applying the second one.

        Args:
            second: The second CliffordTableau to compose with.

        Returns:
            The composed CliffordTableau.

        Raises:
            TypeError: If the type of second is not CliffordTableau.
            ValueError: If the number of qubits in the second tableau mismatch with
                this tableau.
        """
        if not isinstance(second, CliffordTableau):
            raise TypeError("The type for second tableau must be the CliffordTableau type")
        if self.n != second.n:
            raise ValueError(
                f"Mismatched number of qubits of two tableaux: {self.n} vs {second.n}."
            )

        # Convert the underlying data type from bool to int for easier numerical computation.
        m1 = self.matrix().astype(int)
        m2 = second.matrix().astype(int)

        # The following computation is based on Theorem 36 in
        # https://arxiv.org/pdf/2009.03218.pdf.
        # Any pauli string (one stabilizer) in Clifford Tableau should be able to be expressed as
        #    (1i)^p (-1)^s X^(mx) Z^(mz)
        # where p and s are binary scalar and mx and mz are binary vectors.
        num_ys1 = np.sum(m1[:, : self.n] * m1[:, self.n :], axis=1)
        num_ys2 = np.sum(m2[:, : self.n] * m2[:, self.n :], axis=1)

        p1 = np.mod(num_ys1, 2)
        p2 = np.mod(num_ys2, 2)

        # Note the `s` is not equal to `r`, which depends on the number of Y gates.
        # For example, r * Y_1Y_2Y_3 can be expanded into i^3 * r * X_1Z_1 X_2Z_2 X_3Z_3.
        # The global phase is i * (-1) * r ==> s = r + 1 and p = 1.
        s1 = self.rs.astype(int) + np.mod(num_ys1, 4) // 2
        s2 = second.rs.astype(int) + np.mod(num_ys2, 4) // 2

        lmbda = np.zeros((2 * self.n, 2 * self.n))
        lmbda[: self.n, self.n :] = np.eye(self.n)

        m_12 = np.mod(m1 @ m2, 2)
        p_12 = np.mod(p1 + m1 @ p2, 2)
        s_12 = (
            s1
            + m1 @ s2
            + p1 * (m1 @ p2)
            + np.diag(m1 @ np.tril(np.outer(p2, p2.T) + m2 @ lmbda @ m2.T, -1) @ m1.T)
        )
        num_ys12 = np.sum(m_12[:, : self.n] * m_12[:, self.n :], axis=1)
        merged_sign = np.mod(p_12 + 2 * s_12 - num_ys12, 4) // 2

        merged_tableau = CliffordTableau(num_qubits=self.n)
        merged_tableau.xs = m_12[:, : self.n]
        merged_tableau.zs = m_12[:, self.n :]
        merged_tableau.rs = merged_sign

        return merged_tableau
Returns a composed CliffordTableau of this tableau and the second tableau.

        Then composed tableau is equal to (up to global phase) the composed
        unitary operation of the two tableaux, i.e. equivalent to applying the unitary
        operation of this CliffordTableau then applying the second one.

        Args:
            second: The second CliffordTableau to compose with.

        Returns:
            The composed CliffordTableau.

        Raises:
            TypeError: If the type of second is not CliffordTableau.
            ValueError: If the number of qubits in the second tableau mismatch with
                this tableau.
        
cirq-core/cirq/qis/clifford_tableau.py