diff --git a/doc/python_api_reference_vDev.md b/doc/python_api_reference_vDev.md
index 9905cd76c..b922eceab 100644
--- a/doc/python_api_reference_vDev.md
+++ b/doc/python_api_reference_vDev.md
@@ -10366,7 +10366,33 @@ def to_circuit(
                 Note: "graph_state" treats the tableau as a state instead of as a
                 Clifford operation. It will preserve the set of stabilizers, but
                 not the exact choice of generators.
+            "mpp_state": Prepares the tableau's state using MPP and feedback.
+                Gate set: MPP, CX rec, CY rec, CZ rec
+                Circuit qubit count: n
+                Circuit operation count: O(n^2)
+
+                The circuit will be made up of two layers:
+                    1. An MPP layer measuring each of the tableau's stabilizers.
+                    2. A feedback layer using the measurement results to control
+                        whether or not to apply each of the tableau's destabilizers
+                        in order to get the correct sign for each stabilizer.
+
+                Note: "mpp_state" treats the tableau as a state instead of as a
+                Clifford operation. It will preserve the set of stabilizers, but
+                not the exact choice of generators.
+            "mpp_state_unsigned": Prepares the tableau's state up to sign using MPP.
+                Gate set: MPP
+                Circuit qubit count: n
+                Circuit operation count: O(n^2)
 
+                The circuit will contain a series of MPP measurements measuring each
+                of the tableau's stabilizers. The stabilizers are measured in the
+                order used by the tableau (i.e. tableau.z_output(k) is the k'th
+                stabilizer measured).
+
+                Note: "mpp_state_unsigned" treats the tableau as a state instead of
+                as a Clifford operation. It will preserve the set of stabilizers,
+                but not the exact choice of generators.
     Returns:
         The synthesized circuit.
 
@@ -10421,6 +10447,19 @@ def to_circuit(
             H 3
             S 3
         ''')
+
+        >>> tableau.to_circuit("mpp_state_unsigned")
+        stim.Circuit('''
+            MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+        ''')
+
+        >>> tableau.to_circuit("mpp_state")
+        stim.Circuit('''
+            MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+            CX rec[-3] 2 rec[-1] 2
+            CY rec[-4] 0 rec[-3] 0 rec[-3] 3 rec[-2] 3 rec[-1] 0
+            CZ rec[-4] 1 rec[-1] 1
+        ''')
     """
 ```
 
diff --git a/doc/stim.pyi b/doc/stim.pyi
index 6a7ff42f8..32d471208 100644
--- a/doc/stim.pyi
+++ b/doc/stim.pyi
@@ -8070,7 +8070,33 @@ class Tableau:
                     Note: "graph_state" treats the tableau as a state instead of as a
                     Clifford operation. It will preserve the set of stabilizers, but
                     not the exact choice of generators.
+                "mpp_state": Prepares the tableau's state using MPP and feedback.
+                    Gate set: MPP, CX rec, CY rec, CZ rec
+                    Circuit qubit count: n
+                    Circuit operation count: O(n^2)
+
+                    The circuit will be made up of two layers:
+                        1. An MPP layer measuring each of the tableau's stabilizers.
+                        2. A feedback layer using the measurement results to control
+                            whether or not to apply each of the tableau's destabilizers
+                            in order to get the correct sign for each stabilizer.
+
+                    Note: "mpp_state" treats the tableau as a state instead of as a
+                    Clifford operation. It will preserve the set of stabilizers, but
+                    not the exact choice of generators.
+                "mpp_state_unsigned": Prepares the tableau's state up to sign using MPP.
+                    Gate set: MPP
+                    Circuit qubit count: n
+                    Circuit operation count: O(n^2)
 
+                    The circuit will contain a series of MPP measurements measuring each
+                    of the tableau's stabilizers. The stabilizers are measured in the
+                    order used by the tableau (i.e. tableau.z_output(k) is the k'th
+                    stabilizer measured).
+
+                    Note: "mpp_state_unsigned" treats the tableau as a state instead of
+                    as a Clifford operation. It will preserve the set of stabilizers,
+                    but not the exact choice of generators.
         Returns:
             The synthesized circuit.
 
@@ -8125,6 +8151,19 @@ class Tableau:
                 H 3
                 S 3
             ''')
+
+            >>> tableau.to_circuit("mpp_state_unsigned")
+            stim.Circuit('''
+                MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+            ''')
+
+            >>> tableau.to_circuit("mpp_state")
+            stim.Circuit('''
+                MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+                CX rec[-3] 2 rec[-1] 2
+                CY rec[-4] 0 rec[-3] 0 rec[-3] 3 rec[-2] 3 rec[-1] 0
+                CZ rec[-4] 1 rec[-1] 1
+            ''')
         """
     def to_numpy(
         self,
diff --git a/glue/python/src/stim/__init__.pyi b/glue/python/src/stim/__init__.pyi
index 6a7ff42f8..32d471208 100644
--- a/glue/python/src/stim/__init__.pyi
+++ b/glue/python/src/stim/__init__.pyi
@@ -8070,7 +8070,33 @@ class Tableau:
                     Note: "graph_state" treats the tableau as a state instead of as a
                     Clifford operation. It will preserve the set of stabilizers, but
                     not the exact choice of generators.
+                "mpp_state": Prepares the tableau's state using MPP and feedback.
+                    Gate set: MPP, CX rec, CY rec, CZ rec
+                    Circuit qubit count: n
+                    Circuit operation count: O(n^2)
+
+                    The circuit will be made up of two layers:
+                        1. An MPP layer measuring each of the tableau's stabilizers.
+                        2. A feedback layer using the measurement results to control
+                            whether or not to apply each of the tableau's destabilizers
+                            in order to get the correct sign for each stabilizer.
+
+                    Note: "mpp_state" treats the tableau as a state instead of as a
+                    Clifford operation. It will preserve the set of stabilizers, but
+                    not the exact choice of generators.
+                "mpp_state_unsigned": Prepares the tableau's state up to sign using MPP.
+                    Gate set: MPP
+                    Circuit qubit count: n
+                    Circuit operation count: O(n^2)
 
+                    The circuit will contain a series of MPP measurements measuring each
+                    of the tableau's stabilizers. The stabilizers are measured in the
+                    order used by the tableau (i.e. tableau.z_output(k) is the k'th
+                    stabilizer measured).
+
+                    Note: "mpp_state_unsigned" treats the tableau as a state instead of
+                    as a Clifford operation. It will preserve the set of stabilizers,
+                    but not the exact choice of generators.
         Returns:
             The synthesized circuit.
 
@@ -8125,6 +8151,19 @@ class Tableau:
                 H 3
                 S 3
             ''')
+
+            >>> tableau.to_circuit("mpp_state_unsigned")
+            stim.Circuit('''
+                MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+            ''')
+
+            >>> tableau.to_circuit("mpp_state")
+            stim.Circuit('''
+                MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+                CX rec[-3] 2 rec[-1] 2
+                CY rec[-4] 0 rec[-3] 0 rec[-3] 3 rec[-2] 3 rec[-1] 0
+                CZ rec[-4] 1 rec[-1] 1
+            ''')
         """
     def to_numpy(
         self,
diff --git a/src/stim/circuit/gate_target.cc b/src/stim/circuit/gate_target.cc
index 223382183..c2dc17dd4 100644
--- a/src/stim/circuit/gate_target.cc
+++ b/src/stim/circuit/gate_target.cc
@@ -18,6 +18,13 @@
 
 using namespace stim;
 
+GateTarget GateTarget::pauli_xz(uint32_t qubit, bool x, bool z, bool inverted) {
+    if (qubit != (qubit & TARGET_VALUE_MASK)) {
+        throw std::invalid_argument("qubit target larger than " + std::to_string(TARGET_VALUE_MASK));
+    }
+    return {qubit | (TARGET_INVERTED_BIT * inverted) | (TARGET_PAULI_X_BIT * x) | (TARGET_PAULI_Z_BIT * z)};
+}
+
 GateTarget GateTarget::x(uint32_t qubit, bool inverted) {
     if (qubit != (qubit & TARGET_VALUE_MASK)) {
         throw std::invalid_argument("qubit target larger than " + std::to_string(TARGET_VALUE_MASK));
diff --git a/src/stim/circuit/gate_target.h b/src/stim/circuit/gate_target.h
index fa964a386..7a27330fb 100644
--- a/src/stim/circuit/gate_target.h
+++ b/src/stim/circuit/gate_target.h
@@ -39,6 +39,7 @@ struct GateTarget {
     static GateTarget x(uint32_t qubit, bool inverted = false);
     static GateTarget y(uint32_t qubit, bool inverted = false);
     static GateTarget z(uint32_t qubit, bool inverted = false);
+    static GateTarget pauli_xz(uint32_t qubit, bool x, bool z, bool inverted = false);
     static GateTarget qubit(uint32_t qubit, bool inverted = false);
     static GateTarget rec(int32_t lookback);
     static GateTarget sweep_bit(uint32_t index);
diff --git a/src/stim/gates/gates.test.cc b/src/stim/gates/gates.test.cc
index cf5ed6e8a..c2879d2f6 100644
--- a/src/stim/gates/gates.test.cc
+++ b/src/stim/gates/gates.test.cc
@@ -79,8 +79,8 @@ std::pair<std::vector<PauliString<W>>, std::vector<PauliString<W>>> circuit_outp
     }
     TableauSimulator<W> sim1(INDEPENDENT_TEST_RNG(), circuit.count_qubits(), -1);
     TableauSimulator<W> sim2(INDEPENDENT_TEST_RNG(), circuit.count_qubits(), +1);
-    sim1.expand_do_circuit(circuit);
-    sim2.expand_do_circuit(circuit);
+    sim1.safe_do_circuit(circuit);
+    sim2.safe_do_circuit(circuit);
     return {sim1.canonical_stabilizers(), sim2.canonical_stabilizers()};
 }
 
diff --git a/src/stim/simulators/graph_simulator.test.cc b/src/stim/simulators/graph_simulator.test.cc
index eefb74616..3c50ef4e7 100644
--- a/src/stim/simulators/graph_simulator.test.cc
+++ b/src/stim/simulators/graph_simulator.test.cc
@@ -16,8 +16,8 @@ void expect_graph_circuit_is_equivalent(const Circuit &circuit) {
     Circuit converted = sim.to_circuit();
     TableauSimulator<64> sim1(std::mt19937_64{}, n);
     TableauSimulator<64> sim2(std::mt19937_64{}, n);
-    sim1.expand_do_circuit(circuit);
-    sim2.expand_do_circuit(converted);
+    sim1.safe_do_circuit(circuit);
+    sim2.safe_do_circuit(converted);
     auto s1 = sim1.canonical_stabilizers();
     auto s2 = sim2.canonical_stabilizers();
     if (s1 != s2) {
@@ -35,9 +35,9 @@ void expect_graph_sim_effect_matches_tableau_sim(const GraphSimulator &state, co
 
     TableauSimulator<64> tableau_sim1(std::mt19937_64{}, n);
     TableauSimulator<64> tableau_sim2(std::mt19937_64{}, n);
-    tableau_sim1.expand_do_circuit(state.to_circuit());
-    tableau_sim1.expand_do_circuit(effect);
-    tableau_sim2.expand_do_circuit(state_after_effect.to_circuit());
+    tableau_sim1.safe_do_circuit(state.to_circuit());
+    tableau_sim1.safe_do_circuit(effect);
+    tableau_sim2.safe_do_circuit(state_after_effect.to_circuit());
     auto s1 = tableau_sim1.canonical_stabilizers();
     auto s2 = tableau_sim2.canonical_stabilizers();
     if (s1 != s2) {
@@ -240,8 +240,8 @@ TEST(graph_simulator, do_complementation) {
 
         TableauSimulator<64> tableau_sim1(std::mt19937_64{}, 8);
         TableauSimulator<64> tableau_sim2(std::mt19937_64{}, 8);
-        tableau_sim1.expand_do_circuit(sim.to_circuit());
-        tableau_sim2.expand_do_circuit(sim2.to_circuit());
+        tableau_sim1.safe_do_circuit(sim.to_circuit());
+        tableau_sim2.safe_do_circuit(sim2.to_circuit());
         auto s1 = tableau_sim1.canonical_stabilizers();
         auto s2 = tableau_sim2.canonical_stabilizers();
         if (s1 != s2) {
diff --git a/src/stim/simulators/tableau_simulator.h b/src/stim/simulators/tableau_simulator.h
index b67f1d4ab..accab9edb 100644
--- a/src/stim/simulators/tableau_simulator.h
+++ b/src/stim/simulators/tableau_simulator.h
@@ -96,7 +96,7 @@ struct TableauSimulator {
     /// Runs all of the operations in the given circuit.
     ///
     /// Automatically expands the tableau simulator's state, if needed.
-    void expand_do_circuit(const Circuit &circuit, uint64_t reps = 1);
+    void safe_do_circuit(const Circuit &circuit, uint64_t reps = 1);
     void do_operation_ensure_size(const CircuitInstruction &operation);
 
     void apply_tableau(const Tableau<W> &tableau, const std::vector<size_t> &targets);
diff --git a/src/stim/simulators/tableau_simulator.inl b/src/stim/simulators/tableau_simulator.inl
index 4fa7713fc..e3a7f182e 100644
--- a/src/stim/simulators/tableau_simulator.inl
+++ b/src/stim/simulators/tableau_simulator.inl
@@ -1107,7 +1107,7 @@ void TableauSimulator<W>::do_Z(const CircuitInstruction &target_data) {
 template <size_t W>
 simd_bits<W> TableauSimulator<W>::sample_circuit(const Circuit &circuit, std::mt19937_64 &rng, int8_t sign_bias) {
     TableauSimulator<W> sim(std::move(rng), circuit.count_qubits(), sign_bias);
-    sim.expand_do_circuit(circuit);
+    sim.safe_do_circuit(circuit);
 
     const std::vector<bool> &v = sim.measurement_record.storage;
     simd_bits<W> result(v.size());
@@ -1347,7 +1347,7 @@ void TableauSimulator<W>::collapse_isolate_qubit_z(size_t target, TableauTranspo
 }
 
 template <size_t W>
-void TableauSimulator<W>::expand_do_circuit(const Circuit &circuit, uint64_t reps) {
+void TableauSimulator<W>::safe_do_circuit(const Circuit &circuit, uint64_t reps) {
     ensure_large_enough_for_qubits(circuit.count_qubits());
     for (uint64_t k = 0; k < reps; k++) {
         circuit.for_each_operation([&](const CircuitInstruction &op) {
diff --git a/src/stim/simulators/tableau_simulator.pybind.cc b/src/stim/simulators/tableau_simulator.pybind.cc
index f1ce8b978..76f1d835d 100644
--- a/src/stim/simulators/tableau_simulator.pybind.cc
+++ b/src/stim/simulators/tableau_simulator.pybind.cc
@@ -29,10 +29,10 @@ using namespace stim_pybind;
 template <size_t W>
 void do_obj(TableauSimulator<W> &self, const pybind11::object &obj) {
     if (pybind11::isinstance<Circuit>(obj)) {
-        self.expand_do_circuit(pybind11::cast<Circuit>(obj));
+        self.safe_do_circuit(pybind11::cast<Circuit>(obj));
     } else if (pybind11::isinstance<CircuitRepeatBlock>(obj)) {
         const CircuitRepeatBlock &block = pybind11::cast<CircuitRepeatBlock>(obj);
-        self.expand_do_circuit(block.body, block.repeat_count);
+        self.safe_do_circuit(block.body, block.repeat_count);
     } else if (pybind11::isinstance<FlexPauliString>(obj)) {
         const FlexPauliString &pauli_string = pybind11::cast<FlexPauliString>(obj);
         self.ensure_large_enough_for_qubits(pauli_string.value.num_qubits);
@@ -474,7 +474,7 @@ void stim_pybind::pybind_tableau_simulator_methods(
     c.def(
         "do_circuit",
         [](TableauSimulator<MAX_BITWORD_WIDTH> &self, const Circuit &circuit) {
-            self.expand_do_circuit(circuit);
+            self.safe_do_circuit(circuit);
         },
         pybind11::arg("circuit"),
         clean_doc_string(R"DOC(
diff --git a/src/stim/simulators/tableau_simulator.test.cc b/src/stim/simulators/tableau_simulator.test.cc
index 6e4a6bcc9..ffc487997 100644
--- a/src/stim/simulators/tableau_simulator.test.cc
+++ b/src/stim/simulators/tableau_simulator.test.cc
@@ -1425,7 +1425,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, sample_stream_mutates_rng_state, {
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_x, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RX 0
         REPEAT 10000 {
             MX(0.05) 0
@@ -1434,11 +1434,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_x, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MX 0");
+    t.safe_do_circuit("MX 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RX 0 1
         Y 0 1
         REPEAT 5000 {
@@ -1448,13 +1448,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_x, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MX 0");
+    t.safe_do_circuit("MX 0");
     ASSERT_TRUE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_y, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RY 0
         REPEAT 10000 {
             MY(0.05) 0
@@ -1463,11 +1463,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_y, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MY 0");
+    t.safe_do_circuit("MY 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RY 0 1
         X 0 1
         REPEAT 5000 {
@@ -1477,13 +1477,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_y, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MY 0");
+    t.safe_do_circuit("MY 0");
     ASSERT_TRUE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_z, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RZ 0
         REPEAT 10000 {
             MZ(0.05) 0
@@ -1492,11 +1492,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_z, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MZ 0");
+    t.safe_do_circuit("MZ 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RZ 0 1
         X 0 1
         REPEAT 5000 {
@@ -1506,13 +1506,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measurement_z, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MZ 0");
+    t.safe_do_circuit("MZ 0");
     ASSERT_TRUE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_x, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RX 0
         REPEAT 10000 {
             MRX(0.05) 0
@@ -1521,11 +1521,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_x, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MX 0");
+    t.safe_do_circuit("MX 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RX 0 1
         REPEAT 5000 {
             Z 0 1
@@ -1535,13 +1535,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_x, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MX 0");
+    t.safe_do_circuit("MX 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_y, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RY 0 1
         REPEAT 5000 {
             MRY(0.05) 0 1
@@ -1550,11 +1550,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_y, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MY 0");
+    t.safe_do_circuit("MY 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RY 0 1
         REPEAT 5000 {
             X 0 1
@@ -1564,13 +1564,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_y, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MY 0");
+    t.safe_do_circuit("MY 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_z, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RZ 0 1
         REPEAT 5000 {
             MRZ(0.05) 0 1
@@ -1579,11 +1579,11 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_z, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MZ 0");
+    t.safe_do_circuit("MZ 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 
     t.measurement_record.storage.clear();
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         RZ 0 1
         REPEAT 5000 {
             X 0 1
@@ -1593,19 +1593,19 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, noisy_measure_reset_z, {
     m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 10000 - 700);
     ASSERT_LT(m1, 10000 - 300);
-    t.expand_do_circuit("MZ 0");
+    t.safe_do_circuit("MZ 0");
     ASSERT_FALSE(t.measurement_record.storage.back());
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_bad, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    ASSERT_THROW({ t.expand_do_circuit("MPP X0*X0"); }, std::invalid_argument);
-    ASSERT_THROW({ t.expand_do_circuit("MPP X0*Z0"); }, std::invalid_argument);
+    ASSERT_THROW({ t.safe_do_circuit("MPP X0*X0"); }, std::invalid_argument);
+    ASSERT_THROW({ t.safe_do_circuit("MPP X0*Z0"); }, std::invalid_argument);
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_1, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         REPEAT 100 {
             RX 0
             RY 1
@@ -1620,7 +1620,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_1, {
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_4body, {
     for (size_t k = 0; k < 10; k++) {
         TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-        t.expand_do_circuit(R"CIRCUIT(
+        t.safe_do_circuit(R"CIRCUIT(
             MPP X0*X1*X2*X3
             MX 0 1 2 3 4 5
             MPP X2*X3*X4*X5
@@ -1645,7 +1645,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_4body, {
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_epr, {
     for (size_t k = 0; k < 10; k++) {
         TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-        t.expand_do_circuit(R"CIRCUIT(
+        t.safe_do_circuit(R"CIRCUIT(
             MPP X0*X1 Z0*Z1 Y0*Y1
             CNOT 0 1
             H 0
@@ -1665,7 +1665,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_epr, {
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_inversions, {
     for (size_t k = 0; k < 10; k++) {
         TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-        t.expand_do_circuit(R"CIRCUIT(
+        t.safe_do_circuit(R"CIRCUIT(
             MPP !X0*!X1 !X0*X1 X0*!X1 X0*X1 X0 X1 !X0 !X1
         )CIRCUIT");
         auto a = t.measurement_record.storage[0];
@@ -1687,7 +1687,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_inversions, {
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_noisy, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         H 0
         CNOT 0 1
         REPEAT 5000 {
@@ -1697,13 +1697,13 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, measure_pauli_product_noisy, {
     auto m1 = std::accumulate(t.measurement_record.storage.begin(), t.measurement_record.storage.end(), 0);
     ASSERT_GT(m1, 300);
     ASSERT_LT(m1, 700);
-    t.expand_do_circuit("MPP Y0*Y1");
+    t.safe_do_circuit("MPP Y0*Y1");
     ASSERT_EQ(t.measurement_record.storage.back(), true);
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, ignores_sweep_controls, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         X 0
         CNOT sweep[0] 0
         M 0
@@ -1713,7 +1713,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, ignores_sweep_controls, {
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, peek_observable_expectation, {
     TableauSimulator<W> t(INDEPENDENT_TEST_RNG());
-    t.expand_do_circuit(R"CIRCUIT(
+    t.safe_do_circuit(R"CIRCUIT(
         H 0
         CNOT 0 1
         X 0
@@ -2173,15 +2173,15 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, mpad, {
     ASSERT_EQ(sim.inv_state, Tableau<W>(5));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{}));
 
-    sim.expand_do_circuit(Circuit("MPAD 0"));
+    sim.safe_do_circuit(Circuit("MPAD 0"));
     ASSERT_EQ(sim.inv_state, Tableau<W>(5));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{0}));
 
-    sim.expand_do_circuit(Circuit("MPAD 1"));
+    sim.safe_do_circuit(Circuit("MPAD 1"));
     ASSERT_EQ(sim.inv_state, Tableau<W>(5));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{0, 1}));
 
-    sim.expand_do_circuit(Circuit("MPAD 0 0 1 1 0"));
+    sim.safe_do_circuit(Circuit("MPAD 0 0 1 1 0"));
     ASSERT_EQ(sim.inv_state, Tableau<W>(5));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{0, 1, 0, 0, 1, 1, 0}));
 })
@@ -2193,8 +2193,8 @@ void expect_same_final_state(const Tableau<W> &start, const Circuit &c1, const C
     TableauSimulator<W> sim2(INDEPENDENT_TEST_RNG(), n);
     sim1.inv_state = start;
     sim2.inv_state = start;
-    sim1.expand_do_circuit(c1);
-    sim2.expand_do_circuit(c2);
+    sim1.safe_do_circuit(c1);
+    sim2.safe_do_circuit(c2);
     auto t1 = sim1.canonical_stabilizers();
     auto t2 = sim2.canonical_stabilizers();
     if (unsigned_stabilizers) {
@@ -2221,131 +2221,131 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, mxx_myy_mzz_vs_mpp_unsigned, {
 TEST_EACH_WORD_SIZE_W(TableauSimulator, mxx, {
     auto rng = INDEPENDENT_TEST_RNG();
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 5);
-    sim.expand_do_circuit(Circuit("RX 0 1"));
-    sim.expand_do_circuit(Circuit("MXX 0 1"));
+    sim.safe_do_circuit(Circuit("RX 0 1"));
+    sim.safe_do_circuit(Circuit("MXX 0 1"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{false}));
     sim.measurement_record.storage.clear();
 
     sim.inv_state = Tableau<W>::random(5, rng);
-    sim.expand_do_circuit(Circuit("MXX 1 3"));
+    sim.safe_do_circuit(Circuit("MXX 1 3"));
     bool x13 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MXX 1 3"));
-    sim.expand_do_circuit(Circuit("MXX 1 !3"));
-    sim.expand_do_circuit(Circuit("MXX !1 3"));
-    sim.expand_do_circuit(Circuit("MXX !1 !3"));
+    sim.safe_do_circuit(Circuit("MXX 1 3"));
+    sim.safe_do_circuit(Circuit("MXX 1 !3"));
+    sim.safe_do_circuit(Circuit("MXX !1 3"));
+    sim.safe_do_circuit(Circuit("MXX !1 !3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x13, !x13, !x13, x13}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MXX 2 3"));
+    sim.safe_do_circuit(Circuit("MXX 2 3"));
     bool x23 = sim.measurement_record.storage.back();
     bool x12 = x13 ^ x23;
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MXX 1 2"));
+    sim.safe_do_circuit(Circuit("MXX 1 2"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MXX 3 4"));
+    sim.safe_do_circuit(Circuit("MXX 3 4"));
     bool x34 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MXX 1 2 3 4 2 3 1 3"));
+    sim.safe_do_circuit(Circuit("MXX 1 2 3 4 2 3 1 3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12, x34, x23, x13}));
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, myy, {
     auto rng = INDEPENDENT_TEST_RNG();
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 5);
-    sim.expand_do_circuit(Circuit("RY 0 1"));
-    sim.expand_do_circuit(Circuit("MYY 0 1"));
+    sim.safe_do_circuit(Circuit("RY 0 1"));
+    sim.safe_do_circuit(Circuit("MYY 0 1"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{false}));
     sim.measurement_record.storage.clear();
 
     sim.inv_state = Tableau<W>::random(5, rng);
-    sim.expand_do_circuit(Circuit("MYY 1 3"));
+    sim.safe_do_circuit(Circuit("MYY 1 3"));
     bool x13 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MYY 1 3"));
-    sim.expand_do_circuit(Circuit("MYY 1 !3"));
-    sim.expand_do_circuit(Circuit("MYY !1 3"));
-    sim.expand_do_circuit(Circuit("MYY !1 !3"));
+    sim.safe_do_circuit(Circuit("MYY 1 3"));
+    sim.safe_do_circuit(Circuit("MYY 1 !3"));
+    sim.safe_do_circuit(Circuit("MYY !1 3"));
+    sim.safe_do_circuit(Circuit("MYY !1 !3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x13, !x13, !x13, x13}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MYY 2 3"));
+    sim.safe_do_circuit(Circuit("MYY 2 3"));
     bool x23 = sim.measurement_record.storage.back();
     bool x12 = x13 ^ x23;
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MYY 1 2"));
+    sim.safe_do_circuit(Circuit("MYY 1 2"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MYY 3 4"));
+    sim.safe_do_circuit(Circuit("MYY 3 4"));
     bool x34 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MYY 1 2 3 4 2 3 1 3"));
+    sim.safe_do_circuit(Circuit("MYY 1 2 3 4 2 3 1 3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12, x34, x23, x13}));
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, mzz, {
     auto rng = INDEPENDENT_TEST_RNG();
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 5);
-    sim.expand_do_circuit(Circuit("RZ 0 1"));
-    sim.expand_do_circuit(Circuit("MZZ 0 1"));
+    sim.safe_do_circuit(Circuit("RZ 0 1"));
+    sim.safe_do_circuit(Circuit("MZZ 0 1"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{false}));
     sim.measurement_record.storage.clear();
 
     sim.inv_state = Tableau<W>::random(5, rng);
-    sim.expand_do_circuit(Circuit("MZZ 1 3"));
+    sim.safe_do_circuit(Circuit("MZZ 1 3"));
     bool x13 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MZZ 1 3"));
-    sim.expand_do_circuit(Circuit("MZZ 1 !3"));
-    sim.expand_do_circuit(Circuit("MZZ !1 3"));
-    sim.expand_do_circuit(Circuit("MZZ !1 !3"));
+    sim.safe_do_circuit(Circuit("MZZ 1 3"));
+    sim.safe_do_circuit(Circuit("MZZ 1 !3"));
+    sim.safe_do_circuit(Circuit("MZZ !1 3"));
+    sim.safe_do_circuit(Circuit("MZZ !1 !3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x13, !x13, !x13, x13}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MZZ 2 3"));
+    sim.safe_do_circuit(Circuit("MZZ 2 3"));
     bool x23 = sim.measurement_record.storage.back();
     bool x12 = x13 ^ x23;
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MZZ 1 2"));
+    sim.safe_do_circuit(Circuit("MZZ 1 2"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12}));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MZZ 3 4"));
+    sim.safe_do_circuit(Circuit("MZZ 3 4"));
     bool x34 = sim.measurement_record.storage.back();
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit("MZZ 1 2 3 4 2 3 1 3"));
+    sim.safe_do_circuit(Circuit("MZZ 1 2 3 4 2 3 1 3"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{x12, x34, x23, x13}));
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, runs_on_general_circuit, {
     auto circuit = generate_test_circuit_with_all_operations();
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 1);
-    sim.expand_do_circuit(circuit);
+    sim.safe_do_circuit(circuit);
     ASSERT_GT(sim.inv_state.xs.num_qubits, 1);
 })
 
 TEST_EACH_WORD_SIZE_W(TableauSimulator, heralded_erase, {
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 1);
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_ERASE(0) 0 1 2 3 10 11 12 13
     )CIRCUIT"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{0, 0, 0, 0, 0, 0, 0, 0}));
     sim.measurement_record.storage.clear();
     ASSERT_EQ(sim.inv_state, Tableau<W>(14));
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_ERASE(1) 0 1 2 3 4 5 6 10 11 12 13
     )CIRCUIT"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}));
@@ -2356,21 +2356,21 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, heralded_pauli_channel_1, {
     TableauSimulator<W> sim(INDEPENDENT_TEST_RNG(), 1);
     Tableau<W> expected(14);
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_PAULI_CHANNEL_1(0, 0, 0, 0) 0 1 2 3 10 11 12 13
     )CIRCUIT"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{0, 0, 0, 0, 0, 0, 0, 0}));
     ASSERT_EQ(sim.inv_state, Tableau<W>(14));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_PAULI_CHANNEL_1(1, 0, 0, 0) 0 1 2 3 4 5 6 10 11 12 13
     )CIRCUIT"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}));
     ASSERT_EQ(sim.inv_state, Tableau<W>(14));
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_PAULI_CHANNEL_1(0, 1, 0, 0) 13
     )CIRCUIT"));
     ASSERT_EQ(sim.measurement_record.storage, (std::vector<bool>{true}));
@@ -2378,7 +2378,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, heralded_pauli_channel_1, {
     ASSERT_EQ(sim.inv_state, expected);
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_PAULI_CHANNEL_1(0, 0, 1, 0) 5 10
     )CIRCUIT"));
     expected.prepend_Y(5);
@@ -2387,7 +2387,7 @@ TEST_EACH_WORD_SIZE_W(TableauSimulator, heralded_pauli_channel_1, {
     ASSERT_EQ(sim.inv_state, expected);
     sim.measurement_record.storage.clear();
 
-    sim.expand_do_circuit(Circuit(R"CIRCUIT(
+    sim.safe_do_circuit(Circuit(R"CIRCUIT(
         HERALDED_PAULI_CHANNEL_1(0, 0, 0, 1) 1 10 11
     )CIRCUIT"));
     expected.prepend_Z(1);
diff --git a/src/stim/stabilizers/conversions.h b/src/stim/stabilizers/conversions.h
index 4991e2077..a1fb9fc4d 100644
--- a/src/stim/stabilizers/conversions.h
+++ b/src/stim/stabilizers/conversions.h
@@ -133,6 +133,8 @@ Circuit tableau_to_circuit(const Tableau<W> &tableau, const std::string &method)
 template <size_t W>
 Circuit tableau_to_circuit_graph_method(const Tableau<W> &tableau);
 template <size_t W>
+Circuit tableau_to_circuit_mpp_method(const Tableau<W> &tableau, bool skip_sign);
+template <size_t W>
 Circuit tableau_to_circuit_elimination_method(const Tableau<W> &tableau);
 
 /// Converts a unitary matrix into a stabilizer tableau.
diff --git a/src/stim/stabilizers/conversions.inl b/src/stim/stabilizers/conversions.inl
index e93d37428..b6a8edec3 100644
--- a/src/stim/stabilizers/conversions.inl
+++ b/src/stim/stabilizers/conversions.inl
@@ -216,11 +216,17 @@ Circuit tableau_to_circuit(const Tableau<W> &tableau, const std::string &method)
         return tableau_to_circuit_elimination_method(tableau);
     } else if (method == "graph_state") {
         return tableau_to_circuit_graph_method(tableau);
+    } else if (method == "mpp_state") {
+        return tableau_to_circuit_mpp_method(tableau, false);
+    } else if (method == "mpp_state_unsigned") {
+        return tableau_to_circuit_mpp_method(tableau, true);
     } else {
         std::stringstream ss;
         ss << "Unknown method: '" << method << "'. Known methods:\n";
-        ss << "    - 'elimination'";
-        ss << "    - 'graph'";
+        ss << "    - 'elimination'\n";
+        ss << "    - 'graph_state'\n";
+        ss << "    - 'mpp_state'\n";
+        ss << "    - 'mpp_state_unsigned'\n";
         throw std::invalid_argument(ss.str());
     }
 }
@@ -232,6 +238,63 @@ Circuit tableau_to_circuit_graph_method(const Tableau<W> &tableau) {
     return sim.to_circuit(true);
 }
 
+template <size_t W>
+Circuit tableau_to_circuit_mpp_method(const Tableau<W> &tableau, bool skip_sign) {
+    Circuit result;
+    std::vector<GateTarget> targets;
+    size_t n = tableau.num_qubits;
+
+    // Measure each stabilizer with MPP.
+    for (size_t k = 0; k < n; k++) {
+        const auto &stabilizer = tableau.zs[k];
+        bool need_sign = stabilizer.sign;
+        for (size_t q = 0; q < n; q++) {
+            bool x = stabilizer.xs[q];
+            bool z = stabilizer.zs[q];
+            if (x || z) {
+                targets.push_back(GateTarget::pauli_xz(q, x, z, need_sign));
+                targets.push_back(GateTarget::combiner());
+                need_sign = false;
+            }
+        }
+        assert(!targets.empty());
+        targets.pop_back();
+        result.safe_append(GateType::MPP, targets, {});
+        targets.clear();
+    }
+
+    if (!skip_sign) {
+        // Correct each stabilizer's sign with feedback.
+        std::vector<GateTarget> targets_x;
+        std::vector<GateTarget> targets_y;
+        std::vector<GateTarget> targets_z;
+        std::array<std::vector<GateTarget> *, 4> targets_ptrs = {nullptr, &targets_x, &targets_z, &targets_y};
+        for (size_t k = 0; k < n; k++) {
+            const auto &destabilizer = tableau.xs[k];
+            for (size_t q = 0; q < n; q++) {
+                bool x = destabilizer.xs[q];
+                bool z = destabilizer.zs[q];
+                auto *out = targets_ptrs[x + z * 2];
+                if (out != nullptr) {
+                    out->push_back(GateTarget::rec(-(int32_t)(n - k)));
+                    out->push_back(GateTarget::qubit(q));
+                }
+            }
+        }
+        if (!targets_x.empty()) {
+            result.safe_append(GateType::CX, targets_x, {});
+        }
+        if (!targets_y.empty()) {
+            result.safe_append(GateType::CY, targets_y, {});
+        }
+        if (!targets_z.empty()) {
+            result.safe_append(GateType::CZ, targets_z, {});
+        }
+    }
+
+    return result;
+}
+
 template <size_t W>
 Circuit tableau_to_circuit_elimination_method(const Tableau<W> &tableau) {
     Tableau<W> remaining = tableau.inverse();
diff --git a/src/stim/stabilizers/conversions.test.cc b/src/stim/stabilizers/conversions.test.cc
index e14e3c297..7f213e585 100644
--- a/src/stim/stabilizers/conversions.test.cc
+++ b/src/stim/stabilizers/conversions.test.cc
@@ -815,3 +815,59 @@ TEST(conversions, independent_vs_disjoint_xyz_errors_round_trip_fuzz) {
         ASSERT_NEAR(z, z3, 1e-6) << "x=" << x << ",y=" << y << ",z=" << z;
     }
 }
+
+TEST_EACH_WORD_SIZE_W(conversions, fuzz_mpp_circuit_produces_correct_state, {
+    auto rng = INDEPENDENT_TEST_RNG();
+    auto tableau = Tableau<W>::random(10, rng);
+    auto circuit = tableau_to_circuit_mpp_method<W>(tableau, false);
+    TableauSimulator<W> sim(std::move(rng), 10);
+    sim.safe_do_circuit(circuit);
+    auto expected = tableau.stabilizers(true);
+    auto actual = sim.canonical_stabilizers();
+    ASSERT_EQ(actual, expected);
+})
+
+TEST_EACH_WORD_SIZE_W(conversions, fuzz_mpp_circuit_produces_correct_state_unsigned, {
+    auto rng = INDEPENDENT_TEST_RNG();
+    auto tableau = Tableau<W>::random(10, rng);
+    auto circuit = tableau_to_circuit_mpp_method<W>(tableau, true);
+    TableauSimulator<W> sim(std::move(rng), 10);
+    sim.safe_do_circuit(circuit);
+    auto expected = tableau.stabilizers(true);
+    auto actual = sim.canonical_stabilizers();
+    for (auto &e : expected) {
+        e.sign = false;
+    }
+    for (auto &e : actual) {
+        e.sign = false;
+    }
+    ASSERT_EQ(actual, expected);
+})
+
+TEST_EACH_WORD_SIZE_W(conversions, perfect_code_mpp_circuit, {
+    Tableau<W> tableau(5);
+
+    tableau.zs[0] = PauliString<W>("XZZX_");
+    tableau.zs[1] = PauliString<W>("_XZZX");
+    tableau.zs[2] = PauliString<W>("X_XZZ");
+    tableau.zs[3] = PauliString<W>("ZX_XZ");
+    tableau.zs[4] = PauliString<W>("ZZZZZ");
+
+    tableau.xs[0] = PauliString<W>("Z_Z__");
+    tableau.xs[1] = PauliString<W>("ZZZZ_");
+    tableau.xs[2] = PauliString<W>("ZZ_ZZ");
+    tableau.xs[3] = PauliString<W>("_Z__Z");
+    tableau.xs[4] = PauliString<W>("XXXXX");
+
+    ASSERT_TRUE(tableau.satisfies_invariants());
+
+    ASSERT_EQ(tableau_to_circuit_mpp_method<W>(tableau, true), Circuit(R"CIRCUIT(
+        MPP X0*Z1*Z2*X3 X1*Z2*Z3*X4 X0*X2*Z3*Z4 Z0*X1*X3*Z4 Z0*Z1*Z2*Z3*Z4
+    )CIRCUIT"));
+
+    ASSERT_EQ(tableau_to_circuit_mpp_method<W>(tableau, false), Circuit(R"CIRCUIT(
+        MPP X0*Z1*Z2*X3 X1*Z2*Z3*X4 X0*X2*Z3*Z4 Z0*X1*X3*Z4 Z0*Z1*Z2*Z3*Z4
+        CX rec[-1] 0 rec[-1] 1 rec[-1] 2 rec[-1] 3 rec[-1] 4
+        CZ rec[-5] 0 rec[-5] 2 rec[-4] 0 rec[-4] 1 rec[-4] 2 rec[-4] 3 rec[-3] 0 rec[-3] 1 rec[-3] 3 rec[-3] 4 rec[-2] 1 rec[-2] 4
+    )CIRCUIT"));
+})
diff --git a/src/stim/stabilizers/tableau.pybind.cc b/src/stim/stabilizers/tableau.pybind.cc
index 05c15880f..576926e4e 100644
--- a/src/stim/stabilizers/tableau.pybind.cc
+++ b/src/stim/stabilizers/tableau.pybind.cc
@@ -697,7 +697,33 @@ void stim_pybind::pybind_tableau_methods(pybind11::module &m, pybind11::class_<T
                         Note: "graph_state" treats the tableau as a state instead of as a
                         Clifford operation. It will preserve the set of stabilizers, but
                         not the exact choice of generators.
+                    "mpp_state": Prepares the tableau's state using MPP and feedback.
+                        Gate set: MPP, CX rec, CY rec, CZ rec
+                        Circuit qubit count: n
+                        Circuit operation count: O(n^2)
+
+                        The circuit will be made up of two layers:
+                            1. An MPP layer measuring each of the tableau's stabilizers.
+                            2. A feedback layer using the measurement results to control
+                                whether or not to apply each of the tableau's destabilizers
+                                in order to get the correct sign for each stabilizer.
+
+                        Note: "mpp_state" treats the tableau as a state instead of as a
+                        Clifford operation. It will preserve the set of stabilizers, but
+                        not the exact choice of generators.
+                    "mpp_state_unsigned": Prepares the tableau's state up to sign using MPP.
+                        Gate set: MPP
+                        Circuit qubit count: n
+                        Circuit operation count: O(n^2)
 
+                        The circuit will contain a series of MPP measurements measuring each
+                        of the tableau's stabilizers. The stabilizers are measured in the
+                        order used by the tableau (i.e. tableau.z_output(k) is the k'th
+                        stabilizer measured).
+
+                        Note: "mpp_state_unsigned" treats the tableau as a state instead of
+                        as a Clifford operation. It will preserve the set of stabilizers,
+                        but not the exact choice of generators.
             Returns:
                 The synthesized circuit.
 
@@ -752,6 +778,19 @@ void stim_pybind::pybind_tableau_methods(pybind11::module &m, pybind11::class_<T
                     H 3
                     S 3
                 ''')
+
+                >>> tableau.to_circuit("mpp_state_unsigned")
+                stim.Circuit('''
+                    MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+                ''')
+
+                >>> tableau.to_circuit("mpp_state")
+                stim.Circuit('''
+                    MPP X0*Z1*Y2*Y3 !X0*Y1*X2 !Z0*X1*X2*Z3 X0*X1*Z2
+                    CX rec[-3] 2 rec[-1] 2
+                    CY rec[-4] 0 rec[-3] 0 rec[-3] 3 rec[-2] 3 rec[-1] 0
+                    CZ rec[-4] 1 rec[-1] 1
+                ''')
         )DOC")
             .data());
 
diff --git a/src/stim/stabilizers/tableau_pybind_test.py b/src/stim/stabilizers/tableau_pybind_test.py
index 49279d7ac..6a0bc6600 100644
--- a/src/stim/stabilizers/tableau_pybind_test.py
+++ b/src/stim/stabilizers/tableau_pybind_test.py
@@ -953,3 +953,25 @@ def test_to_circuit_graph_state_preserves_stabilizers():
     original = t.to_stabilizers(canonicalize=True)
     reconstructed = c.to_tableau().to_stabilizers(canonicalize=True)
     assert original == reconstructed
+
+
+def test_to_circuit_mpp_preserves_stabilizers():
+    t = stim.Tableau.random(10)
+    original = t.to_stabilizers(canonicalize=True)
+    sim = stim.TableauSimulator()
+    sim.do_circuit(t.to_circuit("mpp_state"))
+    reconstructed = sim.canonical_stabilizers()
+    assert original == reconstructed
+
+
+def test_to_circuit_mpp_unsigned_preserves_stabilizers():
+    t = stim.Tableau.random(10)
+    original = t.to_stabilizers(canonicalize=True)
+    sim = stim.TableauSimulator()
+    sim.do_circuit(t.to_circuit("mpp_state_unsigned"))
+    reconstructed = sim.canonical_stabilizers()
+    for e in original:
+        e.sign = +1
+    for e in reconstructed:
+        e.sign = +1
+    assert original == reconstructed