Add method to get a Clifford gate and phase global phase from unitary

Fixes #2847

Author: maffoo

cirq-core/cirq/ops/clifford_gate.py

@@ -620,6 +620,27 @@ def from_unitary(u: np.ndarray) -> Optional['SingleQubitCliffordGate']:
             _to_clifford_tableau(x_to=x_to, z_to=z_to)
         )
 
+    @classmethod
+    def from_unitary_with_global_phase(cls, u: np.ndarray) -> Optional[Tuple['SingleQubitCliffordGate', complex]]:
+        """Creates Clifford gate with given unitary, including global phase.
+
+        Args:
+            u: 2x2 unitary matrix of a Clifford gate.
+
+        Returns:
+            A tuple of a SingleQubitCliffordGate and a global phase, such that
+            the gate unitary (as given by `cirq.unitary`) times the global phase
+            is identical to the given unitary `u`; or `None` if `u` is not the
+            matrix of a single-qubit Clifford gate.
+        """
+        gate = cls.from_unitary(u)
+        if gate is None:
+            return None
+        # Find the entry with the largest magnitude in the input unitary, to find
+        # the global phase difference between the input unitary and the gate unitary.
+        k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t]))
+        return gate, u[k] / protocols.unitary(gate)[k]
+
     def pauli_tuple(self, pauli: Pauli) -> Tuple[Pauli, bool]:
         """Returns a tuple of a Pauli operator and a boolean.
 

cirq-core/cirq/ops/clifford_gate_test.py

@@ -530,26 +530,38 @@ def test_from_unitary(clifford_gate):
     result_gate = cirq.SingleQubitCliffordGate.from_unitary(u)
     assert result_gate == clifford_gate
 
+    result_gate2, global_phase = cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+    assert result_gate2 == result_gate
+    assert np.allclose(cirq.unitary(result_gate2) * global_phase, u)
+
 
 def test_from_unitary_with_phase_shift():
     u = np.exp(0.42j) * cirq.unitary(cirq.SingleQubitCliffordGate.Y_sqrt)
     gate = cirq.SingleQubitCliffordGate.from_unitary(u)
 
     assert gate == cirq.SingleQubitCliffordGate.Y_sqrt
 
+    gate2, global_phase = cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+    assert gate2 == gate
+    assert np.allclose(cirq.unitary(gate2) * global_phase, u)
+
+
 
 def test_from_unitary_not_clifford():
     # Not a single-qubit gate.
     u = cirq.unitary(cirq.CNOT)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
 
     # Not an unitary matrix.
     u = 2 * cirq.unitary(cirq.X)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
 
     # Not a Clifford gate.
     u = cirq.unitary(cirq.T)
     assert cirq.SingleQubitCliffordGate.from_unitary(u) is None
+    assert cirq.SingleQubitCliffordGate.from_unitary_with_global_phase(u) is None
 
 
 @pytest.mark.parametrize('trans_x,trans_z', _all_rotation_pairs())

cirq-core/cirq/sim/clifford/stabilizer_simulation_state.py

@@ -161,21 +161,15 @@ def _strat_act_from_single_qubit_decompose(
             if not has_unitary(val):
                 return NotImplemented
             u = unitary(val)
-            clifford_gate = SingleQubitCliffordGate.from_unitary(u)
-            if clifford_gate is not None:
-                # Gather the effective unitary applied so as to correct for the
-                # global phase later.
-                final_unitary = np.eye(2)
+            gate_and_phase = SingleQubitCliffordGate.from_unitary_with_global_phase(u)
+            if gate_and_phase is not None:
+                clifford_gate, global_phase = gate_and_phase
+                # Decompose gate into rotations and apply them.
                 for axis, quarter_turns in clifford_gate.decompose_rotation():
                     gate = axis ** (quarter_turns / 2)
                     self._strat_apply_gate(gate, qubits)
-                    final_unitary = np.matmul(unitary(gate), final_unitary)
-
-                # Find the entry with the largest magnitude in the input unitary.
-                k = max(np.ndindex(*u.shape), key=lambda t: abs(u[t]))
-                # Correct the global phase that wasn't conserved in the above
-                # decomposition.
-                self._state.apply_global_phase(u[k] / final_unitary[k])
+                # Apply global phase.
+                self._state.apply_global_phase(global_phase)
                 return True
 
         return NotImplemented