diff --git a/file_lists/perf_files b/file_lists/perf_files
index 5797e069..2d6ac07e 100644
--- a/file_lists/perf_files
+++ b/file_lists/perf_files
@@ -18,4 +18,5 @@ src/stim/stabilizers/tableau.perf.cc
 src/stim/stabilizers/tableau_iter.perf.cc
 src/stim/util_bot/error_decomp.perf.cc
 src/stim/util_bot/probability_util.perf.cc
+src/stim/util_top/reference_sample_tree.perf.cc
 src/stim/util_top/stabilizers_to_tableau.perf.cc
diff --git a/file_lists/source_files_no_main b/file_lists/source_files_no_main
index 0974fe65..8555fcd8 100644
--- a/file_lists/source_files_no_main
+++ b/file_lists/source_files_no_main
@@ -94,5 +94,6 @@ src/stim/util_top/circuit_vs_amplitudes.cc
 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/reference_sample_tree.cc
 src/stim/util_top/simplified_circuit.cc
 src/stim/util_top/transform_without_feedback.cc
diff --git a/file_lists/test_files b/file_lists/test_files
index c3a34b1b..efab5d9b 100644
--- a/file_lists/test_files
+++ b/file_lists/test_files
@@ -89,6 +89,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/reference_sample_tree.test.cc
 src/stim/util_top/simplified_circuit.test.cc
 src/stim/util_top/stabilizers_to_tableau.test.cc
 src/stim/util_top/stabilizers_vs_amplitudes.test.cc
diff --git a/src/stim.h b/src/stim.h
index 4f71b946..d1f55e3a 100644
--- a/src/stim.h
+++ b/src/stim.h
@@ -115,6 +115,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/reference_sample_tree.h"
 #include "stim/util_top/simplified_circuit.h"
 #include "stim/util_top/stabilizers_to_tableau.h"
 #include "stim/util_top/stabilizers_vs_amplitudes.h"
diff --git a/src/stim/io/measure_record.cc b/src/stim/io/measure_record.cc
index 62fcb071..1bdf6cbb 100644
--- a/src/stim/io/measure_record.cc
+++ b/src/stim/io/measure_record.cc
@@ -53,3 +53,22 @@ void MeasureRecord::record_result(bool result) {
     storage.push_back(result);
     unwritten++;
 }
+
+void MeasureRecord::record_results(const std::vector<bool> &results) {
+    storage.insert(storage.end(), results.begin(), results.end());
+    unwritten += results.size();
+}
+
+void MeasureRecord::clear() {
+    unwritten = 0;
+    storage.clear();
+}
+
+void MeasureRecord::discard_results_past_max_lookback() {
+    if (storage.size() > max_lookback) {
+        storage.erase(storage.begin(), storage.end() - max_lookback);
+    }
+    if (unwritten > max_lookback) {
+        unwritten = max_lookback;
+    }
+}
diff --git a/src/stim/io/measure_record.h b/src/stim/io/measure_record.h
index feb7901d..907a3a28 100644
--- a/src/stim/io/measure_record.h
+++ b/src/stim/io/measure_record.h
@@ -47,8 +47,14 @@ struct MeasureRecord {
     /// Args:
     ///     lookback: How far back the measurement is. lookback=1 is the latest measurement, 2 the second latest, etc.
     bool lookback(size_t lookback) const;
+    /// Batch record.
+    void record_results(const std::vector<bool> &results);
     /// Appends a measurement to the record.
     void record_result(bool result);
+    /// Clear the record.
+    void clear();
+    /// Truncates the record to only include bits within the lookback limit.
+    void discard_results_past_max_lookback();
 };
 
 }  // namespace stim
diff --git a/src/stim/simulators/tableau_simulator.perf.cc b/src/stim/simulators/tableau_simulator.perf.cc
index 545e5750..9534be30 100644
--- a/src/stim/simulators/tableau_simulator.perf.cc
+++ b/src/stim/simulators/tableau_simulator.perf.cc
@@ -14,6 +14,7 @@
 
 #include "stim/simulators/tableau_simulator.h"
 
+#include "stim/gen/circuit_gen_params.h"
 #include "stim/perf.perf.h"
 
 using namespace stim;
diff --git a/src/stim/util_top/reference_sample_tree.cc b/src/stim/util_top/reference_sample_tree.cc
new file mode 100644
index 00000000..bb8764aa
--- /dev/null
+++ b/src/stim/util_top/reference_sample_tree.cc
@@ -0,0 +1,207 @@
+#include "stim/util_top/reference_sample_tree.h"
+
+using namespace stim;
+
+bool ReferenceSampleTree::empty() const {
+    if (repetitions == 0) {
+        return true;
+    }
+    if (!prefix_bits.empty()) {
+        return false;
+    }
+    for (const auto &child: suffix_children) {
+        if (!child.empty()) {
+            return false;
+        }
+    }
+    return true;
+}
+
+void ReferenceSampleTree::flatten_and_simplify_into(std::vector<ReferenceSampleTree> &out) const {
+    if (repetitions == 0) {
+        return;
+    }
+
+    // Flatten children.
+    std::vector<ReferenceSampleTree> flattened;
+    if (!prefix_bits.empty()) {
+        flattened.push_back(ReferenceSampleTree{
+            .prefix_bits = prefix_bits,
+            .suffix_children = {},
+            .repetitions = 1,
+        });
+    }
+    for (const auto &child : suffix_children) {
+        child.flatten_and_simplify_into(flattened);
+    }
+
+    // Fuse children.
+    std::vector<ReferenceSampleTree> fused;
+    if (!flattened.empty()) {
+        fused.push_back(std::move(flattened[0]));
+    }
+    for (size_t k = 1; k < flattened.size(); k++) {
+        auto &dst = fused.back();
+        auto &src = flattened[k];
+
+        // Combine children with identical contents by adding their rep counts.
+        if (dst.prefix_bits == src.prefix_bits && dst.suffix_children == src.suffix_children) {
+            dst.repetitions += src.repetitions;
+
+        // Fuse children with unrepeated contents if they can be fused.
+        } else if (src.repetitions == 1 && dst.repetitions == 1 && dst.suffix_children.empty()) {
+            dst.suffix_children = std::move(src.suffix_children);
+            dst.prefix_bits.insert(dst.prefix_bits.end(), src.prefix_bits.begin(), src.prefix_bits.end());
+
+        } else {
+            fused.push_back(std::move(src));
+        }
+    }
+
+    if (repetitions == 1) {
+        // Un-nest all the children.
+        for (auto &e : fused) {
+            out.push_back(e);
+        }
+    } else if (fused.size() == 1) {
+        // Merge with single child.
+        fused[0].repetitions *= repetitions;
+        out.push_back(std::move(fused[0]));
+    } else if (fused.empty()) {
+        // Nothing to report.
+    } else if (fused[0].suffix_children.empty() && fused[0].repetitions == 1) {
+        // Take payload from first child.
+        ReferenceSampleTree result = std::move(fused[0]);
+        fused.erase(fused.begin());
+        result.repetitions = repetitions;
+        result.suffix_children = std::move(fused);
+        out.push_back(std::move(result));
+    } else {
+        out.push_back(ReferenceSampleTree{
+            .prefix_bits={},
+            .suffix_children=std::move(fused),
+            .repetitions=repetitions,
+        });
+    }
+}
+
+/// Finds how far back feedback operations ever look, within the loop.
+uint64_t stim::max_feedback_lookback_in_loop(const Circuit &loop) {
+    uint64_t furthest_lookback = 0;
+    for (const auto &inst : loop.operations) {
+        if (inst.gate_type == GateType::REPEAT) {
+            furthest_lookback = std::max(furthest_lookback, max_feedback_lookback_in_loop(inst.repeat_block_body(loop)));
+        } else {
+            auto f = GATE_DATA[inst.gate_type].flags;
+            if ((f & GateFlags::GATE_CAN_TARGET_BITS) && (f & GateFlags::GATE_TARGETS_PAIRS)) {
+                // Feedback-capable operation. Check for any measurement record targets.
+                for (auto t : inst.targets) {
+                    if (t.is_measurement_record_target()) {
+                        furthest_lookback = std::max(furthest_lookback, (uint64_t)-t.rec_offset());
+                    }
+                }
+            }
+        }
+    }
+    return furthest_lookback;
+}
+
+void ReferenceSampleTree::try_factorize(size_t period_factor) {
+    if (prefix_bits.size() != 0 || suffix_children.size() % period_factor != 0) {
+        return;
+    }
+
+    // Check if contents are periodic with the factor.
+    size_t h = suffix_children.size() / period_factor;
+    for (size_t k = h; k < suffix_children.size(); k++) {
+        if (suffix_children[k - h] != suffix_children[k]) {
+            return;
+        }
+    }
+
+    // Factorize.
+    suffix_children.resize(h);
+    repetitions *= period_factor;
+}
+
+ReferenceSampleTree ReferenceSampleTree::simplified() const {
+    std::vector<ReferenceSampleTree> flat;
+    flatten_and_simplify_into(flat);
+    if (flat.empty()) {
+        return ReferenceSampleTree();
+    } else if (flat.size() == 1) {
+        return std::move(flat[0]);
+    }
+
+    ReferenceSampleTree result;
+    result.repetitions = 1;
+
+     // Take payload from first child.
+    if (flat[0].repetitions == 1 && flat[0].suffix_children.empty()) {
+         result = std::move(flat[0]);
+         flat.erase(flat.begin());
+    }
+
+    result.suffix_children = std::move(flat);
+    return result;
+}
+
+size_t ReferenceSampleTree::size() const {
+    size_t result = prefix_bits.size();
+    for (const auto &child: suffix_children) {
+        result += child.size();
+    }
+    return result * repetitions;
+}
+
+void ReferenceSampleTree::decompress_into(std::vector<bool> &output) const {
+    for (uint64_t k = 0; k < repetitions; k++) {
+        output.insert(output.end(), prefix_bits.begin(), prefix_bits.end());
+        for (const auto &child: suffix_children) {
+            child.decompress_into(output);
+        }
+    }
+}
+
+ReferenceSampleTree ReferenceSampleTree::from_circuit_reference_sample(const Circuit &circuit) {
+    auto stats = circuit.compute_stats();
+    std::mt19937_64 irrelevant_rng{0};
+    CompressedReferenceSampleHelper<MAX_BITWORD_WIDTH> helper(
+        TableauSimulator<MAX_BITWORD_WIDTH>(
+            std::move(irrelevant_rng),
+            stats.num_qubits,
+            +1,
+            MeasureRecord(stats.max_lookback)));
+    return helper.do_loop_with_tortoise_hare_folding(circuit, 1).simplified();
+}
+
+std::string ReferenceSampleTree::str() const {
+    std::stringstream ss;
+    ss << *this;
+    return ss.str();
+}
+
+bool ReferenceSampleTree::operator==(const ReferenceSampleTree &other) const {
+    return repetitions == other.repetitions &&
+           prefix_bits == other.prefix_bits &&
+           suffix_children == other.suffix_children;
+}
+bool ReferenceSampleTree::operator!=(const ReferenceSampleTree &other) const {
+    return !(*this == other);
+}
+
+std::ostream &stim::operator<<(std::ostream &out, const ReferenceSampleTree &v) {
+    out << v.repetitions << "*";
+    out << "(";
+    out << "'";
+    for (auto b : v.prefix_bits) {
+        out << "01"[b];
+    }
+    out << "'";
+    for (const auto &child : v.suffix_children) {
+        out << "+";
+        out << child;
+    }
+    out << ")";
+    return out;
+}
diff --git a/src/stim/util_top/reference_sample_tree.h b/src/stim/util_top/reference_sample_tree.h
new file mode 100644
index 00000000..2b439b61
--- /dev/null
+++ b/src/stim/util_top/reference_sample_tree.h
@@ -0,0 +1,86 @@
+#ifndef _STIM_UTIL_TOP_REFERENCE_SAMPLE_TREE_H
+#define _STIM_UTIL_TOP_REFERENCE_SAMPLE_TREE_H
+
+#include "stim/simulators/tableau_simulator.h"
+
+namespace stim {
+
+/// A compressed tree representation of a reference sample.
+struct ReferenceSampleTree {
+    /// Raw bits to output before bits from the children.
+    std::vector<bool> prefix_bits;
+    /// Compressed representations of additional bits to output after the prefix.
+    std::vector<ReferenceSampleTree> suffix_children;
+    /// The number of times to repeatedly output the prefix+suffix bits.
+    size_t repetitions = 0;
+
+    /// Initializes a reference sample tree containing a reference sample for the given circuit.
+    static ReferenceSampleTree from_circuit_reference_sample(const Circuit &circuit);
+
+    /// Returns a tree with the same compressed contents, but a simpler tree structure.
+    ReferenceSampleTree simplified() const;
+
+    /// Checks if two trees are exactly the same, including structure (not just uncompressed contents).
+    bool operator==(const ReferenceSampleTree &other) const;
+    /// Checks if two trees are not exactly the same, including structure (not just uncompressed contents).
+    bool operator!=(const ReferenceSampleTree &other) const;
+    /// Returns a simple description of the tree's structure, like "5*('101'+6*('11'))".
+    std::string str() const;
+
+    /// Determines whether the tree contains any bits at all.
+    bool empty() const;
+    /// Computes the total size of the uncompressed bits represented by the tree.
+    size_t size() const;
+
+    /// Writes the contents of the tree into the given output vector.
+    void decompress_into(std::vector<bool> &output) const;
+
+    /// Folds redundant children into the repetition count, if they repeat this many times.
+    ///
+    /// For example, if the tree's children are [A, B, C, A, B, C] and the tree has no
+    /// prefix, then `try_factorize(2)` will reduce the children to [A, B, C] and double
+    /// the repetition count.
+    void try_factorize(size_t period_factor);
+
+   private:
+    /// Helper method for `simplified`.
+    void flatten_and_simplify_into(std::vector<ReferenceSampleTree> &out) const;
+};
+std::ostream &operator<<(std::ostream &out, const ReferenceSampleTree &v);
+
+/// Helper class for computing compressed reference samples.
+template <size_t W>
+struct CompressedReferenceSampleHelper {
+    TableauSimulator<W> sim;
+
+    CompressedReferenceSampleHelper(TableauSimulator<MAX_BITWORD_WIDTH> sim) : sim(sim) {
+    }
+
+    /// Processes a loop with no top-level folding.
+    ///
+    /// Loops containing within the body of this loop (or circuit body) may
+    /// still be compressed. Only the top-level loop is not folded.
+    ReferenceSampleTree do_loop_with_no_folding(
+        const Circuit &loop,
+        uint64_t reps);
+
+    /// Runs tortoise-and-hare analysis of the loop while simulating its
+    /// reference sample, in order to attempt to return a compressed
+    /// representation.
+    ReferenceSampleTree do_loop_with_tortoise_hare_folding(
+        const Circuit &loop,
+        uint64_t reps);
+
+    bool in_same_recent_state_as(
+        const CompressedReferenceSampleHelper<W> &other,
+        uint64_t max_record_lookback,
+        bool allow_false_negative) const;
+};
+
+uint64_t max_feedback_lookback_in_loop(const Circuit &loop);
+
+}  // namespace stim
+
+#include "stim/util_top/reference_sample_tree.inl"
+
+#endif
diff --git a/src/stim/util_top/reference_sample_tree.inl b/src/stim/util_top/reference_sample_tree.inl
new file mode 100644
index 00000000..746920dc
--- /dev/null
+++ b/src/stim/util_top/reference_sample_tree.inl
@@ -0,0 +1,142 @@
+#include "stim/util_top/reference_sample_tree.h"
+
+namespace stim {
+
+template <size_t W>
+ReferenceSampleTree CompressedReferenceSampleHelper<W>::do_loop_with_no_folding(const Circuit &loop, uint64_t reps) {
+    ReferenceSampleTree result;
+    result.repetitions = 1;
+    size_t start_size = sim.measurement_record.storage.size();
+
+    auto flush_recorded_into_result = [&]() {
+        size_t end_size = sim.measurement_record.storage.size();
+        if (end_size > start_size) {
+            result.suffix_children.push_back({});
+            auto &child = result.suffix_children.back();
+            child.repetitions = 1;
+            child.prefix_bits.insert(
+                child.prefix_bits.end(),
+                sim.measurement_record.storage.begin() + start_size,
+                sim.measurement_record.storage.begin() + end_size);
+        }
+        start_size = end_size;
+    };
+
+    for (size_t k = 0; k < reps; k++) {
+        for (const auto &inst : loop.operations) {
+            if (inst.gate_type == GateType::REPEAT) {
+                uint64_t repeats = inst.repeat_block_rep_count();
+                const auto& block = inst.repeat_block_body(loop);
+                flush_recorded_into_result();
+                result.suffix_children.push_back(do_loop_with_tortoise_hare_folding(block, repeats));
+                start_size = sim.measurement_record.storage.size();
+            } else {
+                sim.do_gate(inst);
+            }
+        }
+    }
+
+    flush_recorded_into_result();
+    return result;
+}
+
+template <size_t W>
+ReferenceSampleTree CompressedReferenceSampleHelper<W>::do_loop_with_tortoise_hare_folding(const Circuit &loop, uint64_t reps) {
+    if (reps < 10) {
+        // Probably not worth the overhead of tortoise-and-hare. Just run it raw.
+        return do_loop_with_no_folding(loop, reps);
+    }
+
+    ReferenceSampleTree result;
+    result.repetitions = 1;
+
+    CompressedReferenceSampleHelper<W> tortoise(sim);
+    CompressedReferenceSampleHelper<W> hare(std::move(sim));
+    uint64_t max_feedback_lookback = max_feedback_lookback_in_loop(loop);
+    uint64_t tortoise_steps = 0;
+    uint64_t hare_steps = 0;
+    while (hare_steps < reps) {
+        hare_steps++;
+        result.suffix_children.push_back(hare.do_loop_with_no_folding(loop, 1));
+        assert(result.suffix_children.size() == hare_steps);
+
+        if (tortoise.in_same_recent_state_as(hare, max_feedback_lookback, hare_steps < 10)) {
+            break;
+        }
+
+        // Tortoise advances half as quickly.
+        if (hare_steps & 1) {
+            tortoise_steps++;
+            tortoise.do_loop_with_no_folding(loop, 1);
+        }
+    }
+
+    if (hare_steps == reps) {
+        // No periodic state found before reaching the end of the loop.
+        sim = std::move(hare.sim);
+        return result;
+    }
+
+    // Run more loop iterations until the remaining iterations are a multiple of the found period.
+    assert(result.suffix_children.size() == hare_steps);
+    uint64_t period = hare_steps - tortoise_steps;
+    size_t period_steps_left = (reps - hare_steps) / period;
+    while ((reps - hare_steps) % period) {
+        result.suffix_children.push_back(hare.do_loop_with_no_folding(loop, 1));
+        hare_steps += 1;
+    }
+    assert(hare_steps + period_steps_left * period == reps);
+    assert(hare_steps >= period);
+    sim = std::move(hare.sim);
+
+    // Move the periodic measurements out of the hare's tail, into a loop node.
+    ReferenceSampleTree loop_contents;
+    for (size_t k = hare_steps - period; k < hare_steps; k++) {
+        loop_contents.suffix_children.push_back(std::move(result.suffix_children[k]));
+    }
+    result.suffix_children.resize(hare_steps - period);
+
+    // Add skipped iterations' measurement data into the sim's measurement record.
+    loop_contents.repetitions = 1;
+    sim.measurement_record.discard_results_past_max_lookback();
+    for (size_t k = 0; k < period_steps_left && sim.measurement_record.storage.size() < sim.measurement_record.max_lookback * 2; k++) {
+        loop_contents.decompress_into(sim.measurement_record.storage);
+    }
+    sim.measurement_record.discard_results_past_max_lookback();
+
+    // Add the loop node to the output data.
+    loop_contents.repetitions = period_steps_left + 1;
+    loop_contents.try_factorize(2);
+    loop_contents.try_factorize(3);
+    loop_contents.try_factorize(5);
+    result.suffix_children.push_back(std::move(loop_contents));
+
+    return result;
+}
+
+template <size_t W>
+bool CompressedReferenceSampleHelper<W>::in_same_recent_state_as(
+    const CompressedReferenceSampleHelper<W> &other,
+    uint64_t max_record_lookback,
+    bool allow_false_negative) const {
+    const auto &s1 = sim.measurement_record.storage;
+    const auto &s2 = other.sim.measurement_record.storage;
+
+    // Check that recent measurements gave identical results.
+    if (s1.size() < max_record_lookback || s2.size() < max_record_lookback) {
+        return false;
+    }
+    for (size_t k = 0; k < max_record_lookback; k++) {
+        if (s1[s1.size() - k - 1] != s2[s2.size() - k - 1]) {
+            return false;
+        }
+    }
+
+    // Check that quantum states are identical.
+    if (allow_false_negative) {
+        return sim.inv_state == other.sim.inv_state;
+    }
+    return sim.canonical_stabilizers() == other.sim.canonical_stabilizers();
+}
+
+}  // namespace stim
diff --git a/src/stim/util_top/reference_sample_tree.perf.cc b/src/stim/util_top/reference_sample_tree.perf.cc
new file mode 100644
index 00000000..c245f157
--- /dev/null
+++ b/src/stim/util_top/reference_sample_tree.perf.cc
@@ -0,0 +1,52 @@
+#include "stim/util_top/reference_sample_tree.h"
+
+#include "stim/gen/circuit_gen_params.h"
+#include "stim/gen/gen_surface_code.h"
+#include "stim/perf.perf.h"
+
+using namespace stim;
+
+BENCHMARK(reference_sample_tree_surface_code_d31_r1000000000) {
+    CircuitGenParameters params(1000000000, 31, "rotated_memory_x");
+    auto circuit = generate_surface_code_circuit(params).circuit;
+    simd_bits<MAX_BITWORD_WIDTH> ref(0);
+    auto total = 0;
+    benchmark_go([&]() {
+        auto result = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+        total += result.empty();
+    })
+        .goal_millis(25);
+    if (total) {
+        std::cerr << "data dependence";
+    }
+}
+
+BENCHMARK(reference_sample_tree_nested_circuit) {
+    Circuit circuit(R"CIRCUIT(
+        M 0
+        REPEAT 100000 {
+            REPEAT 100000 {
+                REPEAT 100000 {
+                    X 0
+                    M 0
+                }
+                X 0
+                M 0
+            }
+            X 0
+            M 0
+        }
+        X 0
+        M 0
+    )CIRCUIT");
+    simd_bits<MAX_BITWORD_WIDTH> ref(0);
+    auto total = 0;
+    benchmark_go([&]() {
+        auto result = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+        total += result.empty();
+    })
+        .goal_micros(230);
+    if (total) {
+        std::cerr << "data dependence";
+    }
+}
diff --git a/src/stim/util_top/reference_sample_tree.test.cc b/src/stim/util_top/reference_sample_tree.test.cc
new file mode 100644
index 00000000..858a8e66
--- /dev/null
+++ b/src/stim/util_top/reference_sample_tree.test.cc
@@ -0,0 +1,268 @@
+#include "stim/util_top/reference_sample_tree.h"
+
+#include "gtest/gtest.h"
+
+#include "stim/gen/gen_surface_code.h"
+
+using namespace stim;
+
+void expect_tree_matches_normal_reference_sample_of(const ReferenceSampleTree &tree, const Circuit &circuit) {
+    std::vector<bool> decompressed;
+    tree.decompress_into(decompressed);
+    simd_bits<MAX_BITWORD_WIDTH> actual(decompressed.size());
+    for (size_t k = 0; k < decompressed.size(); k++) {
+        actual[k] = decompressed[k];
+    }
+    auto expected = TableauSimulator<MAX_BITWORD_WIDTH>::reference_sample_circuit(circuit);
+    EXPECT_EQ(actual, expected);
+}
+
+TEST(ReferenceSampleTree, equality) {
+    ReferenceSampleTree empty1{
+        .prefix_bits={},
+        .suffix_children={},
+        .repetitions=0,
+    };
+    ReferenceSampleTree empty2;
+    ASSERT_EQ(empty1, empty2);
+
+    ASSERT_FALSE(empty1 != empty2);
+    ASSERT_NE(empty1, (ReferenceSampleTree{.prefix_bits={},.suffix_children{},.repetitions=1}));
+    ASSERT_NE(empty1, (ReferenceSampleTree{.prefix_bits={0},.suffix_children{},.repetitions=0}));
+    ASSERT_NE(empty1, (ReferenceSampleTree{.prefix_bits={},.suffix_children{{}},.repetitions=0}));
+}
+
+TEST(ReferenceSampleTree, str) {
+    ASSERT_EQ((ReferenceSampleTree{
+        .prefix_bits={},
+        .suffix_children={},
+        .repetitions=0,
+    }.str()), "0*('')");
+
+    ASSERT_EQ((ReferenceSampleTree{
+        .prefix_bits={1, 1, 0, 1},
+        .suffix_children={},
+        .repetitions=0,
+    }.str()), "0*('1101')");
+
+    ASSERT_EQ((ReferenceSampleTree{
+        .prefix_bits={1, 1, 0, 1},
+        .suffix_children={},
+        .repetitions=2,
+    }.str()), "2*('1101')");
+
+    ASSERT_EQ((ReferenceSampleTree{
+        .prefix_bits={1, 1, 0, 1},
+        .suffix_children={
+            ReferenceSampleTree{
+                .prefix_bits={1},
+                .suffix_children={},
+                .repetitions=5,
+            }
+        },
+        .repetitions=2,
+    }.str()), "2*('1101'+5*('1'))");
+}
+
+TEST(ReferenceSampleTree, simplified) {
+    ReferenceSampleTree raw{
+        .prefix_bits={},
+        .suffix_children={
+            ReferenceSampleTree{
+                .prefix_bits={},
+                .suffix_children={},
+                .repetitions=1,
+            },
+            ReferenceSampleTree{
+                .prefix_bits={1, 0, 1},
+                .suffix_children={{}},
+                .repetitions=0,
+            },
+            ReferenceSampleTree{
+                .prefix_bits={1, 1, 1},
+                .suffix_children={},
+                .repetitions=2,
+            },
+        },
+        .repetitions=3,
+    };
+    ASSERT_EQ(raw.simplified().str(), "6*('111')");
+}
+
+TEST(ReferenceSampleTree, decompress_into) {
+    std::vector<bool> result;
+    ReferenceSampleTree{
+        .prefix_bits={1, 1, 0, 1},
+        .suffix_children={
+            ReferenceSampleTree{
+                .prefix_bits={1},
+                .suffix_children={},
+                .repetitions=5,
+            }
+        },
+        .repetitions=2,
+    }.decompress_into(result);
+    ASSERT_EQ(result, (std::vector<bool>{1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1}));
+
+    result.clear();
+    ReferenceSampleTree{
+        .prefix_bits={1, 1, 0, 1},
+        .suffix_children={
+            ReferenceSampleTree{
+                .prefix_bits={1, 0, 1},
+                .suffix_children={},
+                .repetitions=8,
+            },
+            ReferenceSampleTree{
+                .prefix_bits={0, 0},
+                .suffix_children={},
+                .repetitions=1,
+            },
+        },
+        .repetitions=1,
+    }.decompress_into(result);
+    ASSERT_EQ(result, (std::vector<bool>{1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0}));
+}
+
+TEST(ReferenceSampleTree, simple_circuit) {
+    Circuit circuit(R"CIRCUIT(
+        M 0
+        X 0
+        M 0
+    )CIRCUIT");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+    ASSERT_EQ(ref.str(), "1*('01')");
+}
+
+TEST(ReferenceSampleTree, simple_loop) {
+    Circuit circuit(R"CIRCUIT(
+        REPEAT 50 {
+            M 0
+            X 0
+            M 0
+        }
+    )CIRCUIT");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+    ASSERT_EQ(ref.str(), "25*('0110')");
+}
+
+TEST(ReferenceSampleTree, period4_loop) {
+    Circuit circuit(R"CIRCUIT(
+        M 0
+        X 0
+        M 0
+        REPEAT 50 {
+            CX 0 1 1 2 2 3
+            M 0 1 2 3
+        }
+        X 0
+        M 0
+        X 2
+        M 2 2 2 2 2
+        MPAD 1 0 1 0 1 1
+    )CIRCUIT");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    ASSERT_EQ(ref.size(), circuit.count_measurements());
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+    ASSERT_EQ(ref.str(), "1*('01111110101100100011111010'+11*('1100100011111010')+1*('000000101011'))");
+}
+
+TEST(ReferenceSampleTree, feedback) {
+    Circuit circuit(R"CIRCUIT(
+        MPAD 0 0 1 0
+        REPEAT 200 {
+            CX rec[-4] 1
+            M 1
+        }
+    )CIRCUIT");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    ASSERT_EQ(ref.size(), circuit.count_measurements());
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+    ASSERT_EQ(ref.str(), "1*('0010'+2*('0')+4*('1')+1*('01011001000111')+12*('101011001000111'))");
+}
+
+TEST(max_feedback_lookback_in_loop, simple) {
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit()), 0);
+
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit(R"CIRCUIT(
+        REPEAT 100 {
+            REPEAT 100 {
+                M 0
+                X 0
+                M 0
+            }
+            REPEAT 200 {
+                M 0
+                DETECTOR rec[-1]
+            }
+            X 1
+            CX 1 0
+        }
+    )CIRCUIT")), 0);
+
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit(R"CIRCUIT(
+        CX rec[-1] 0
+    )CIRCUIT")), 1);
+
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit(R"CIRCUIT(
+        CZ 0 rec[-2]
+    )CIRCUIT")), 2);
+
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit(R"CIRCUIT(
+        CZ 0 rec[-2]
+        CY 0 rec[-3]
+    )CIRCUIT")), 3);
+
+    ASSERT_EQ(max_feedback_lookback_in_loop(Circuit(R"CIRCUIT(
+        CZ 0 rec[-2]
+        REPEAT 100 {
+            CX rec[-5] 0
+        }
+    )CIRCUIT")), 5);
+}
+
+TEST(ReferenceSampleTree, nested_loops) {
+    Circuit circuit(R"CIRCUIT(
+        REPEAT 100 {
+            REPEAT 100 {
+                M 0
+                X 0
+                M 0
+            }
+            REPEAT 200 {
+                M 0
+            }
+            X 1
+            CX 1 0
+        }
+    )CIRCUIT");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+    ASSERT_EQ(ref.str(), "1*(''+50*('0110')+200*('0')+50*('1001')+200*('1')+50*('1001')+200*('1')+50*('0110')+200*('0')+24*(''+50*('0110')+200*('0')+50*('1001')+200*('1')+50*('1001')+200*('1')+50*('0110')+200*('0')))");
+}
+
+TEST(ReferenceSampleTree, surface_code) {
+    CircuitGenParameters params(10000, 5, "rotated_memory_x");
+    auto circuit = generate_surface_code_circuit(params).circuit;
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    ASSERT_EQ(ref.str(), "1*(''+10000*('000000000000000000000000')+1*('0000000000000000000000000'))");
+}
+
+TEST(ReferenceSampleTree, surface_code_with_pauli) {
+    CircuitGenParameters params(10000, 3, "rotated_memory_x");
+    auto circuit = generate_surface_code_circuit(params).circuit;
+    circuit.blocks[0].append_from_text("X 10 11 12 13");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    ASSERT_EQ(ref.str(), "1*(''+2*('00000000')+4999*('0110000000110000')+1*('000000000'))");
+}
+
+TEST(ReferenceSampleTree, surface_code_with_pauli_vs_normal_reference_sample) {
+    CircuitGenParameters params(20, 3, "rotated_memory_x");
+    auto circuit = generate_surface_code_circuit(params).circuit;
+    circuit.blocks[0].append_from_text("X 10 11 12 13");
+    auto ref = ReferenceSampleTree::from_circuit_reference_sample(circuit);
+    ASSERT_EQ(ref.size(), circuit.count_measurements());
+    expect_tree_matches_normal_reference_sample_of(ref, circuit);
+}