Fix the CZ ops on the Clifford Tableau Implementation (#4182)

Implementation idea of CZ is `CZ(q0, q1) = H(q1)CNOT(q0, q1)H(q1)` so the operations are just composition of these three.
The problem happens at `rs` of tableau. Original implementation seemed merge the `rs` of `H` and `CNOT` together. But it is not equivalent. See the fix code.

FYI, I discovered this because I keep encountering the failure in my Clifford decomposition test. I didn't realize it is the CZ implementation in the cirq has problem for a long time :(

Author: BichengYing

cirq-core/cirq/ops/common_gates.py

@@ -1079,8 +1079,11 @@ def _act_on_(self, args: 'cirq.ActOnArgs', qubits: Sequence['cirq.Qid']):
                     tableau.zs[:, q2].copy(),
                     tableau.xs[:, q2].copy(),
                 )
+                tableau.rs[:] ^= tableau.xs[:, q2] & tableau.zs[:, q2]
                 tableau.rs[:] ^= (
-                    tableau.xs[:, q1] & tableau.zs[:, q2] & (tableau.xs[:, q2] ^ tableau.zs[:, q1])
+                    tableau.xs[:, q1]
+                    & tableau.zs[:, q2]
+                    & (~(tableau.xs[:, q2] ^ tableau.zs[:, q1]))
                 )
                 tableau.xs[:, q2] ^= tableau.xs[:, q1]
                 tableau.zs[:, q1] ^= tableau.zs[:, q2]

cirq-core/cirq/ops/common_gates_test.py

@@ -532,6 +532,32 @@ def test_cz_act_on_tableau():
         cirq.act_on(cirq.CZ ** 1.5, args, cirq.LineQubit.range(2))
 
 
+def test_cz_act_on_equivalent_to_h_cx_h_tableau():
+    args1 = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=2),
+        qubits=cirq.LineQubit.range(2),
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    args2 = cirq.ActOnCliffordTableauArgs(
+        tableau=cirq.CliffordTableau(num_qubits=2),
+        qubits=cirq.LineQubit.range(2),
+        prng=np.random.RandomState(),
+        log_of_measurement_results={},
+    )
+    cirq.act_on(cirq.S, args=args1, qubits=[cirq.LineQubit(1)], allow_decompose=False)
+    cirq.act_on(cirq.S, args=args2, qubits=[cirq.LineQubit(1)], allow_decompose=False)
+
+    # Args1 uses H*CNOT*H
+    cirq.act_on(cirq.H, args=args1, qubits=[cirq.LineQubit(1)], allow_decompose=False)
+    cirq.act_on(cirq.CNOT, args=args1, qubits=cirq.LineQubit.range(2), allow_decompose=False)
+    cirq.act_on(cirq.H, args=args1, qubits=[cirq.LineQubit(1)], allow_decompose=False)
+    # Args2 uses CZ
+    cirq.act_on(cirq.CZ, args=args2, qubits=cirq.LineQubit.range(2), allow_decompose=False)
+
+    assert args1.tableau == args2.tableau
+
+
 foo = sympy.Symbol('foo')