Add MBQC decompositions to the gates documentation (#942)

- Fix `stim.Flow.__init__` not xor-sorting measurements/observables
- Add `stim.Flow.__mul__`

Author: Strilanc

doc/gates.md

@@ -220,6 +220,7 @@ API references for stable versions are kept on the [stim github wiki](https://gi
 - [`stim.Flow`](#stim.Flow)
     - [`stim.Flow.__eq__`](#stim.Flow.__eq__)
     - [`stim.Flow.__init__`](#stim.Flow.__init__)
+    - [`stim.Flow.__mul__`](#stim.Flow.__mul__)
     - [`stim.Flow.__ne__`](#stim.Flow.__ne__)
     - [`stim.Flow.__repr__`](#stim.Flow.__repr__)
     - [`stim.Flow.__str__`](#stim.Flow.__str__)
@@ -9165,6 +9166,40 @@ def __init__(
     """
 ```
 
+<a name="stim.Flow.__mul__"></a>
+```python
+# stim.Flow.__mul__
+
+# (in class stim.Flow)
+def __mul__(
+    self,
+    rhs: stim.Flow,
+) -> stim.Flow:
+    """Computes the product of two flows.
+
+    Args:
+        rhs: The right hand side of the multiplication.
+
+    Returns:
+        The product of the two flows.
+
+    Raises:
+        ValueError: The inputs anti-commute (their product would be anti-Hermitian).
+            For example, 1 -> X times 1 -> Y fails because it would give 1 -> iZ.
+
+    Examples:
+        >>> import stim
+        >>> stim.Flow("X -> X") * stim.Flow("Z -> Z")
+        stim.Flow("Y -> Y")
+
+        >>> stim.Flow("1 -> XX") * stim.Flow("1 -> ZZ")
+        stim.Flow("1 -> -YY")
+
+        >>> stim.Flow("X -> rec[-1]") * stim.Flow("X -> rec[-2]")
+        stim.Flow("_ -> rec[-2] xor rec[-1]")
+    """
+```
+
 <a name="stim.Flow.__ne__"></a>
 ```python
 # stim.Flow.__ne__

doc/python_api_reference_vDev.md

@@ -7309,6 +7309,33 @@ class Flow:
             >>> stim.Flow("X -> Y xor obs[3] xor obs[3] xor obs[3]")
             stim.Flow("X -> Y xor obs[3]")
         """
+    def __mul__(
+        self,
+        rhs: stim.Flow,
+    ) -> stim.Flow:
+        """Computes the product of two flows.
+
+        Args:
+            rhs: The right hand side of the multiplication.
+
+        Returns:
+            The product of the two flows.
+
+        Raises:
+            ValueError: The inputs anti-commute (their product would be anti-Hermitian).
+                For example, 1 -> X times 1 -> Y fails because it would give 1 -> iZ.
+
+        Examples:
+            >>> import stim
+            >>> stim.Flow("X -> X") * stim.Flow("Z -> Z")
+            stim.Flow("Y -> Y")
+
+            >>> stim.Flow("1 -> XX") * stim.Flow("1 -> ZZ")
+            stim.Flow("1 -> -YY")
+
+            >>> stim.Flow("X -> rec[-1]") * stim.Flow("X -> rec[-2]")
+            stim.Flow("_ -> rec[-2] xor rec[-1]")
+        """
     def __ne__(
         self,
         arg0: stim.Flow,

doc/stim.pyi

@@ -95,6 +95,7 @@ src/stim/util_top/export_crumble_url.cc
 src/stim/util_top/export_qasm.cc
 src/stim/util_top/export_quirk_url.cc
 src/stim/util_top/has_flow.cc
+src/stim/util_top/mbqc_decomposition.cc
 src/stim/util_top/reference_sample_tree.cc
 src/stim/util_top/simplified_circuit.cc
 src/stim/util_top/transform_without_feedback.cc

file_lists/source_files_no_main

@@ -92,6 +92,7 @@ src/stim/util_top/export_crumble_url.test.cc
 src/stim/util_top/export_qasm.test.cc
 src/stim/util_top/export_quirk_url.test.cc
 src/stim/util_top/has_flow.test.cc
+src/stim/util_top/mbqc_decomposition.test.cc
 src/stim/util_top/reference_sample_tree.test.cc
 src/stim/util_top/simplified_circuit.test.cc
 src/stim/util_top/stabilizers_to_tableau.test.cc

file_lists/test_files

@@ -7309,6 +7309,33 @@ class Flow:
             >>> stim.Flow("X -> Y xor obs[3] xor obs[3] xor obs[3]")
             stim.Flow("X -> Y xor obs[3]")
         """
+    def __mul__(
+        self,
+        rhs: stim.Flow,
+    ) -> stim.Flow:
+        """Computes the product of two flows.
+
+        Args:
+            rhs: The right hand side of the multiplication.
+
+        Returns:
+            The product of the two flows.
+
+        Raises:
+            ValueError: The inputs anti-commute (their product would be anti-Hermitian).
+                For example, 1 -> X times 1 -> Y fails because it would give 1 -> iZ.
+
+        Examples:
+            >>> import stim
+            >>> stim.Flow("X -> X") * stim.Flow("Z -> Z")
+            stim.Flow("Y -> Y")
+
+            >>> stim.Flow("1 -> XX") * stim.Flow("1 -> ZZ")
+            stim.Flow("1 -> -YY")
+
+            >>> stim.Flow("X -> rec[-1]") * stim.Flow("X -> rec[-2]")
+            stim.Flow("_ -> rec[-2] xor rec[-1]")
+        """
     def __ne__(
         self,
         arg0: stim.Flow,

glue/python/src/stim/__init__.pyi

@@ -117,6 +117,7 @@
 #include "stim/util_top/export_qasm.h"
 #include "stim/util_top/export_quirk_url.h"
 #include "stim/util_top/has_flow.h"
+#include "stim/util_top/mbqc_decomposition.h"
 #include "stim/util_top/reference_sample_tree.h"
 #include "stim/util_top/simplified_circuit.h"
 #include "stim/util_top/stabilizers_to_tableau.h"

src/stim.h

@@ -19,23 +19,24 @@
 #include <iostream>
 #include <map>
 #include <set>
-#include <stim/circuit/circuit.h>
-
-#include "command_convert.h"
-#include "command_detect.h"
-#include "command_diagram.h"
-#include "command_explain_errors.h"
-#include "command_gen.h"
-#include "command_m2d.h"
-#include "command_repl.h"
-#include "command_sample.h"
-#include "command_sample_dem.h"
+
+#include "stim/circuit/circuit.h"
+#include "stim/cmd/command_convert.h"
+#include "stim/cmd/command_detect.h"
+#include "stim/cmd/command_diagram.h"
+#include "stim/cmd/command_explain_errors.h"
+#include "stim/cmd/command_gen.h"
+#include "stim/cmd/command_m2d.h"
+#include "stim/cmd/command_repl.h"
+#include "stim/cmd/command_sample.h"
+#include "stim/cmd/command_sample_dem.h"
 #include "stim/cmd/command_analyze_errors.h"
 #include "stim/gates/gates.h"
 #include "stim/io/stim_data_formats.h"
 #include "stim/stabilizers/flow.h"
 #include "stim/stabilizers/tableau.h"
 #include "stim/util_bot/arg_parse.h"
+#include "stim/util_top/mbqc_decomposition.h"
 
 using namespace stim;
 
@@ -255,13 +256,39 @@ void print_decomposition(Acc &out, const Gate &gate) {
         out << "# The following circuit is equivalent (up to global phase) to `";
         out << undecomposed.str() << "`";
         out << decomposition;
-        if (Circuit(decomposition) == Circuit(undecomposed.str())) {
+        Circuit c(decomposition);
+        if (c == Circuit(undecomposed.str())) {
             out << "\n# (The decomposition is trivial because this gate is in the target gate set.)\n";
+        } else if (c.operations.empty()) {
+            out << "\n# (The decomposition is empty because this gate has no effect.)\n";
         }
         out.change_indent(-4);
     }
 }
 
+void print_mbqc_decomposition(Acc &out, const Gate &gate) {
+    const char *decomposition = mbqc_decomposition(gate.id);
+    if (decomposition != nullptr) {
+        std::stringstream undecomposed;
+        auto decomp_targets = gate_decomposition_help_targets_for_gate_type(gate.id);
+        undecomposed << CircuitInstruction{gate.id, {}, decomp_targets, ""};
+
+        out << "MBQC Decomposition (into MX, MY, MZ, MXX, MZZ, and Pauli feedback):\n";
+        out.change_indent(+4);
+        out << "# The following circuit performs `";
+        out << undecomposed.str() << "` (but affects the measurement record and an ancilla qubit)";
+        out << decomposition;
+        Circuit c(decomposition);
+        if (c == Circuit(undecomposed.str())) {
+            out << "\n# (The decomposition is trivial because this gate is in the target gate set.)\n";
+        } else if (c.operations.empty()) {
+            out << "\n# (The decomposition is empty because this gate has no effect.)\n";
+        }
+
+        out.change_indent(-4);
+    }
+}
+
 void print_stabilizer_generators(Acc &out, const Gate &gate) {
     auto flows = gate.flows<MAX_BITWORD_WIDTH>();
     if (flows.empty()) {
@@ -425,6 +452,7 @@ std::string generate_per_gate_help_markdown(const Gate &alt_gate, int indent, bo
     print_bloch_vector(out, gate);
     print_unitary_matrix(out, gate);
     print_decomposition(out, gate);
+    print_mbqc_decomposition(out, gate);
     out.flush();
     return out.settled;
 }

src/stim/cmd/command_help.cc

@@ -36,7 +36,11 @@ struct Flow {
     /// ENTRIES MUST BE SORTED AND UNIQUE.
     std::vector<uint32_t> observables;
 
+    /// Fixes non-unique non-sorted measurements and observables.
+    void canonicalize();
+
     static Flow<W> from_str(std::string_view text);
+    Flow<W> operator*(const Flow<W> &rhs) const;
     bool operator<(const Flow<W> &other) const;
     bool operator==(const Flow<W> &other) const;
     bool operator!=(const Flow<W> &other) const;

src/stim/stabilizers/flow.h

@@ -258,4 +258,34 @@ std::ostream &operator<<(std::ostream &out, const Flow<W> &flow) {
     return out;
 }
 
+template <size_t W>
+void Flow<W>::canonicalize() {
+    size_t measurements_kept = inplace_xor_sort(SpanRef<int32_t>(measurements)).size();
+    size_t observables_kept = inplace_xor_sort(SpanRef<uint32_t>(observables)).size();
+    measurements.resize(measurements_kept);
+    observables.resize(observables_kept);
+}
+
+template <size_t W>
+Flow<W> Flow<W>::operator*(const Flow<W> &rhs) const {
+    Flow<W> result = *this;
+
+    result.input.ensure_num_qubits(rhs.input.num_qubits, 1.1);
+    result.output.ensure_num_qubits(rhs.output.num_qubits, 1.1);
+    uint8_t log_i = 0;
+    log_i -= result.input.ref().inplace_right_mul_returning_log_i_scalar(rhs.input.ref());
+    log_i += result.output.ref().inplace_right_mul_returning_log_i_scalar(rhs.output.ref());
+    if (log_i & 1) {
+        throw std::invalid_argument(str() + " anticommutes with " + rhs.str());
+    }
+    if (log_i & 2) {
+        result.output.sign ^= true;
+    }
+
+    result.measurements.insert(result.measurements.end(), rhs.measurements.begin(), rhs.measurements.end());
+    result.observables.insert(result.observables.end(), rhs.observables.begin(), rhs.observables.end());
+    result.canonicalize();
+    return result;
+}
+
 }  // namespace stim

src/stim/stabilizers/flow.inl

@@ -101,14 +101,13 @@ static Flow<MAX_BITWORD_WIDTH> py_init_flow(
                     result.measurements.push_back(pybind11::cast<int32_t>(h));
                 }
             }
-            inplace_xor_sort(SpanRef<int32_t>(result.measurements));
         }
         if (!included_observables.is_none()) {
             for (const auto &h : included_observables) {
                 result.observables.push_back(pybind11::cast<uint32_t>(h));
             }
-            inplace_xor_sort(SpanRef<uint32_t>(result.observables));
         }
+        result.canonicalize();
         return result;
     }
 
@@ -287,6 +286,36 @@ void stim_pybind::pybind_flow_methods(pybind11::module &m, pybind11::class_<Flow
         )DOC")
             .data());
 
+    c.def(
+        "__mul__",
+        &Flow<MAX_BITWORD_WIDTH>::operator*,
+        pybind11::arg("rhs"),
+        clean_doc_string(R"DOC(
+            Computes the product of two flows.
+
+            Args:
+                rhs: The right hand side of the multiplication.
+
+            Returns:
+                The product of the two flows.
+
+            Raises:
+                ValueError: The inputs anti-commute (their product would be anti-Hermitian).
+                    For example, 1 -> X times 1 -> Y fails because it would give 1 -> iZ.
+
+            Examples:
+                >>> import stim
+                >>> stim.Flow("X -> X") * stim.Flow("Z -> Z")
+                stim.Flow("Y -> Y")
+
+                >>> stim.Flow("1 -> XX") * stim.Flow("1 -> ZZ")
+                stim.Flow("1 -> -YY")
+
+                >>> stim.Flow("X -> rec[-1]") * stim.Flow("X -> rec[-2]")
+                stim.Flow("_ -> rec[-2] xor rec[-1]")
+        )DOC")
+            .data());
+
     c.def(pybind11::self == pybind11::self, "Determines if two flows have identical contents.");
     c.def(pybind11::self != pybind11::self, "Determines if two flows have non-identical contents.");
     c.def("__str__", &Flow<MAX_BITWORD_WIDTH>::str, "Returns a shorthand description of the flow.");

src/stim/stabilizers/flow.pybind.cc

@@ -176,3 +176,24 @@ def test_obs_include_pauli_terms_sensitivity():
     assert zs == 0
     assert 256 <= ys <= 768
     assert xs == ys
+
+
+def test_flow_canonicalization():
+    assert stim.Flow(measurements=[4, 0, 4]) == stim.Flow(measurements=[0])
+    assert stim.Flow(included_observables=[4, 0, 4]) == stim.Flow(included_observables=[0])
+
+
+def test_flow_multiplication():
+    assert stim.Flow("XYZ -> 1") * stim.Flow("1 -> XYZ") == stim.Flow("XYZ -> XYZ")
+    assert stim.Flow("XX_ -> 1") * stim.Flow("_XX -> 1") == stim.Flow("X_X -> 1")
+    assert stim.Flow("1 -> XX_") * stim.Flow("1 -> _XX") == stim.Flow("1 -> X_X")
+    assert stim.Flow("1 -> rec[-1] xor rec[-3]") * stim.Flow("1 -> rec[-1] xor rec[-2]") == stim.Flow("1 -> rec[-2] xor rec[-3]")
+    assert stim.Flow("1 -> obs[1] xor obs[3]") * stim.Flow("1 -> obs[1] xor obs[2]") == stim.Flow("1 -> obs[2] xor obs[3]")
+    assert stim.Flow("X -> X") * stim.Flow("Z -> Z") == stim.Flow("Y -> Y")
+    assert stim.Flow("1 -> XX") * stim.Flow("1 -> ZZ") == stim.Flow("1 -> -YY")
+    assert stim.Flow("1 -> obs[1]") * stim.Flow("1 -> obs[1]") == stim.Flow("1 -> 1")
+    assert stim.Flow("1 -> rec[1]") * stim.Flow("1 -> rec[1]") == stim.Flow("1 -> 1")
+    with pytest.raises(ValueError, match="anticommute"):
+        _ = stim.Flow("1 -> X") * stim.Flow("1 -> Y")
+    with pytest.raises(ValueError, match="anticommute"):
+        _ = stim.Flow("1 -> Y") * stim.Flow("1 -> X")

src/stim/stabilizers/flow_pybind_test.py

@@ -114,14 +114,14 @@ PauliStringRef<W> &PauliStringRef<W>::operator*=(const PauliStringRef<W> &rhs) {
 
 template <size_t W>
 uint8_t PauliStringRef<W>::inplace_right_mul_returning_log_i_scalar(const PauliStringRef<W> &rhs) noexcept {
-    assert(num_qubits == rhs.num_qubits);
+    assert(num_qubits >= rhs.num_qubits);
 
     // Accumulator registers for counting mod 4 in parallel across each bit position.
     simd_word<W> cnt1{};
     simd_word<W> cnt2{};
 
-    xs.for_each_word(
-        zs, rhs.xs, rhs.zs, [&cnt1, &cnt2](simd_word<W> &x1, simd_word<W> &z1, simd_word<W> &x2, simd_word<W> &z2) {
+    rhs.xs.for_each_word(
+        rhs.zs, xs, zs, [&cnt1, &cnt2](simd_word<W> &x2, simd_word<W> &z2, simd_word<W> &x1, simd_word<W> &z1) {
             // Update the left hand side Paulis.
             auto old_x1 = x1;
             auto old_z1 = z1;