Add CPP gate

doc/gates.md

@@ -78,6 +78,7 @@
     - [MYY](#MYY)
     - [MZZ](#MZZ)
 - Generalized Pauli Product Gates
+    - [CPP](#CPP)
     - [MPP](#MPP)
     - [SPP](#SPP)
     - [SPP_DAG](#SPP_DAG)
@@ -2976,6 +2977,72 @@ Decomposition (into H, S, CX, M, R):
 
 ## Generalized Pauli Product Gates
 
+<a name="CPP"></a>
+### The 'CPP' Instruction
+
+The generalized CNOT gate. Negates states in the intersection of Pauli product observables.
+
+Parens Arguments:
+
+    This instruction takes no parens arguments.
+
+Targets:
+
+    A series of pairs of Pauli products to intersect.
+
+    Each Pauli product is a series of Pauli targets (`[XYZ]#`), record targets (`rec[-#]`),
+    or sweep targets (`sweep[#]`) separated by combiners (`*`). Each product can be negated
+    by prefixing a Pauli target in the product with an inverter (`!`).
+
+    The number of products must be even. CPP X1 Y2 Z3 isn't allowed.
+    Within each pair of products, the pair must commute. CPP X1 Z1 isn't allowed.
+
+Examples:
+
+    # Perform a CNOT gate with qubit 1 as the control and qubit 2 as the target.
+    CPP X1 Z2
+
+    # Perform a CZ gate between qubit 2 and qubit 5, then between qubit 3 and 4.
+    CPP Z2 Z5 Z3 Z4
+
+    # Perform many CX gates, all controlled by qubit 2, targeting qubits 5 through 10.
+    CPP Z2 X5*X6*X7*X8*X9*X10
+
+    # Swap qubit 1 and qubit 5 by negating their overlap with the singlet state.
+    CPP X1*X5 Z1*Z5
+
+    # Negate the amplitude of the |00> state.
+    CPP !Z0 !Z1
+
+    # Measure qubit 0 and do Pauli operations conditioned on the measurement returning TRUE.
+    M 0
+    CPP rec[-1] X1*Y2*Z3
+
+Stabilizer Generators (for `CPP X0*Y1 Z2*Z3`):
+
+    X___ -> X___
+    Z___ -> Z_ZZ
+    _X__ -> _XZZ
+    _Z__ -> _ZZZ
+    __X_ -> XYX_
+    __Z_ -> __Z_
+    ___X -> XY_X
+    ___Z -> ___Z
+
+Decomposition (into H, S, CX, M, R):
+
+    # The following circuit is equivalent (up to global phase) to `CPP X0*Y1 Z2*Z3`
+    CX 3 2
+    CX 1 0
+    S 1
+    S 1
+    S 1
+    CX 2 1
+    S 1
+    CX 1 0
+    CX 3 2
+
+
 <a name="MPP"></a>
 ### The 'MPP' Instruction
 

glue/crumble/test/generated_gate_name_list.test.js

@@ -70,6 +70,7 @@ RZ
 MXX
 MYY
 MZZ
+CPP
 MPP
 SPP
 SPP_DAG

src/stim/circuit/circuit.test.cc

@@ -207,6 +207,33 @@ TEST(circuit, parse_mpp) {
     ASSERT_EQ(c.operations[0].targets.size(), 3);
 }
 
+TEST(circuit, parse_cpp) {
+    ASSERT_THROW({ Circuit("CPP X1"); }, std::invalid_argument);
+    ASSERT_THROW({ Circuit("CPP 1 2"); }, std::invalid_argument);
+    ASSERT_THROW({ Circuit("CPP X1 X2 X3"); }, std::invalid_argument);
+    ASSERT_THROW({ Circuit("CPP X1*X2 X2 X3"); }, std::invalid_argument);
+    ASSERT_THROW({ Circuit("CPP X1*X4 X2*Z5*Z9*Z10 X3"); }, std::invalid_argument);
+    ASSERT_THROW({ Circuit("CPP rec[-1]"); }, std::invalid_argument);
+
+    Circuit c;
+
+    c = Circuit("CPP");
+    ASSERT_EQ(c.operations.size(), 1);
+
+    c = Circuit("CPP rec[-1] X1*Y2*Z3");
+    ASSERT_EQ(c.operations.size(), 1);
+
+    c = Circuit("CPP sweep[0] rec[-1]*X1*sweep[2]");
+    ASSERT_EQ(c.operations.size(), 1);
+
+    c = Circuit("CPP X1 Z2");
+    ASSERT_EQ(c.operations.size(), 1);
+    ASSERT_EQ(c.operations[0].targets.size(), 2);
+    ASSERT_EQ(
+        c.operations[0].targets,
+        ((SpanRef<const GateTarget>)std::vector<GateTarget>{GateTarget::x(1), GateTarget::z(2)}));
+}
+
 TEST(circuit, parse_spp) {
     ASSERT_THROW({ Circuit("SPP 1"); }, std::invalid_argument);
     ASSERT_THROW({ Circuit("SPP rec[-1]"); }, std::invalid_argument);
@@ -1705,6 +1732,7 @@ Circuit stim::generate_test_circuit_with_all_operations() {
 
         # Pauli Product Gates
         MPP X0*Y1*Z2 Z0*Z1
+        CPP X3*X4*X5 Z3*Z4*Y6 Y7 Y8
         SPP X0*Y1*Z2 X3
         SPP_DAG X0*Y1*Z2 X2
         TICK

src/stim/circuit/export_qasm.cc

@@ -138,6 +138,14 @@ struct QasmExporter {
         out << "// --- end decomposed MPP\n";
     }
 
+    void output_decomposed_cpp_operation(const CircuitInstruction &inst) {
+        out << "// --- begin decomposed " << inst << "\n";
+        decompose_cpp_operation_with_reverse_independence(inst, stats.num_qubits, [&](const CircuitInstruction &inst) {
+            output_instruction(inst);
+        });
+        out << "// --- end decomposed CPP\n";
+    }
+
     void output_decomposed_spp_or_spp_dag_operation(const CircuitInstruction &inst) {
         out << "// --- begin decomposed " << inst << "\n";
         decompose_spp_or_spp_dag_operation(inst, stats.num_qubits, false, [&](const CircuitInstruction &inst) {
@@ -528,6 +536,9 @@ struct QasmExporter {
             case GateType::SPP_DAG:
                 output_decomposed_spp_or_spp_dag_operation(instruction);
                 return;
+            case GateType::CPP:
+                output_decomposed_cpp_operation(instruction);
+                return;
 
             default:
                 break;

src/stim/circuit/export_qasm.test.cc

@@ -403,6 +403,30 @@ cx q[1], q[0];
 measure q[0] -> rec[3];
 cx q[1], q[0];
 // --- end decomposed MPP
+// --- begin decomposed CPP X3*X4*X5 Z3*Z4*Y6 Y7 Y8
+h q[3];
+h q[4];
+h q[5];
+cx q[4], q[3];
+cx q[5], q[3];
+h q[4];
+hyz q[6];
+cx q[6], q[4];
+cz q[3], q[4];
+cx q[6], q[4];
+hyz q[6];
+h q[4];
+cx q[5], q[3];
+cx q[4], q[3];
+h q[3];
+h q[4];
+h q[5];
+hyz q[7];
+hyz q[8];
+cz q[7], q[8];
+hyz q[7];
+hyz q[8];
+// --- end decomposed CPP
 // --- begin decomposed SPP X0*Y1*Z2 X3
 h q[0];
 hyz q[1];
@@ -611,6 +635,24 @@ cx q[1], q[0];
 measure q[0] -> rec[3];
 cx q[1], q[0];
 // --- end decomposed MPP
+hyz q[6];
+cx q[6], q[3];
+cx q[4], q[3];
+h q[4];
+h q[5];
+cx q[5], q[4];
+cz q[3], q[4];
+cx q[5], q[4];
+h q[4];
+h q[5];
+cx q[4], q[3];
+cx q[6], q[3];
+hyz q[8];
+hyz q[7];
+hyz q[6];
+cz q[8], q[7];
+hyz q[8];
+hyz q[7];
 // --- begin decomposed SPP X0*Y1*Z2 X3
 h q[0];
 hyz q[1];

src/stim/circuit/gate_decomposition.cc

@@ -269,6 +269,156 @@ void stim::decompose_spp_or_spp_dag_operation(
     }
 }
 
+static void decompose_cpp_operation_with_reverse_independence_helper(
+    CircuitInstruction cpp_op,
+    PauliStringRef<64> obs1,
+    PauliStringRef<64> obs2,
+    std::span<const GateTarget> classical_bits1,
+    std::span<const GateTarget> classical_bits2,
+    const std::function<void(const CircuitInstruction &inst)> &do_instruction_callback,
+    Circuit *workspace,
+    std::vector<GateTarget> *buf) {
+    assert(obs1.num_qubits == obs2.num_qubits);
+
+    if (!obs1.commutes(obs2)) {
+        std::stringstream ss;
+        ss << "Attempted to CPP two anticommuting observables.\n";
+        ss << "    obs1: " << obs1 << "\n";
+        ss << "    obs2: " << obs2 << "\n";
+        ss << "    instruction: " << cpp_op;
+        throw std::invalid_argument(ss.str());
+    }
+
+    workspace->clear();
+    auto apply_fixup = [&](CircuitInstruction inst) {
+        workspace->safe_append(inst);
+        obs1.do_instruction(inst);
+        obs2.do_instruction(inst);
+    };
+
+    auto reduce = [&](PauliStringRef<64> target_obs) {
+        // Turn all non-identity terms into Z terms.
+        for_each_active_qubit_in<true, false>(target_obs, [&](uint32_t q) {
+            GateTarget t = GateTarget::qubit(q);
+            apply_fixup({target_obs.zs[q] ? GateType::H_YZ : GateType::H, {}, &t});
+        });
+
+        // Cancel any extra Z terms.
+        uint64_t pivot = UINT64_MAX;
+        for_each_active_qubit_in<true, true>(target_obs, [&](uint32_t q) {
+            if (pivot == UINT64_MAX) {
+                pivot = q;
+            } else {
+                std::array<GateTarget, 2> ts{GateTarget::qubit(q), GateTarget::qubit(pivot)};
+                apply_fixup({GateType::CX, {}, ts});
+            }
+        });
+
+        return pivot;
+    };
+
+    uint64_t pivot1 = reduce(obs1);
+    uint64_t pivot2 = reduce(obs2);
+
+    if (pivot1 == pivot2 && pivot1 != UINT64_MAX) {
+        // Both observables had identical quantum parts (up to sign).
+        // If their sign differed, we should do nothing.
+        // If their sign matched, we should apply Z to obs1.
+        assert(obs1.xs == obs2.xs);
+        assert(obs1.zs == obs2.zs);
+        obs2.zs[pivot2] = false;
+        obs2.sign ^= obs1.sign;
+        obs2.sign ^= true;
+        pivot2 = UINT64_MAX;
+    }
+    assert(obs1.weight() <= 1);
+    assert(obs2.weight() <= 1);
+    assert((pivot1 == UINT64_MAX) == (obs1.weight() == 0));
+    assert((pivot2 == UINT64_MAX) == (obs2.weight() == 0));
+    assert(pivot1 == UINT64_MAX || obs1.xs[pivot1] + 2 * obs1.zs[pivot1] == 2);
+    assert(pivot1 == UINT64_MAX || obs2.xs[pivot1] + 2 * obs2.zs[pivot1] == 0);
+    assert(pivot2 == UINT64_MAX || obs1.xs[pivot2] + 2 * obs1.zs[pivot2] == 0);
+    assert(pivot2 == UINT64_MAX || obs2.xs[pivot2] + 2 * obs2.zs[pivot2] == 2);
+
+    // Apply rewrites.
+    workspace->for_each_operation(do_instruction_callback);
+
+    // Handle the quantum-quantum interaction.
+    if (pivot1 != UINT64_MAX && pivot2 != UINT64_MAX) {
+        assert(pivot1 != pivot2);
+        std::array<GateTarget, 2> ts{GateTarget::qubit(pivot1), GateTarget::qubit(pivot2)};
+        do_instruction_callback({GateType::CZ, {}, ts});
+    }
+
+    // Handle sign and classical feedback into obs1.
+    if (pivot1 != UINT64_MAX) {
+        for (const auto &t : classical_bits2) {
+            std::array<GateTarget, 2> ts{t, GateTarget::qubit(pivot1)};
+            do_instruction_callback({GateType::CZ, {}, ts});
+        }
+        if (obs2.sign) {
+            GateTarget t = GateTarget::qubit(pivot1);
+            do_instruction_callback({GateType::Z, {}, &t});
+        }
+    }
+
+    // Handle sign and classical feedback into obs2.
+    if (pivot2 != UINT64_MAX) {
+        for (const auto &t : classical_bits1) {
+            std::array<GateTarget, 2> ts{t, GateTarget::qubit(pivot2)};
+            do_instruction_callback({GateType::CZ, {}, ts});
+        }
+        if (obs1.sign) {
+            GateTarget t = GateTarget::qubit(pivot2);
+            do_instruction_callback({GateType::Z, {}, &t});
+        }
+    }
+
+    // Undo rewrites.
+    workspace->for_each_operation_reverse([&](CircuitInstruction inst) {
+        assert(inst.args.empty());
+        if (inst.gate_type == GateType::CX) {
+            buf->clear();
+            for (size_t k = inst.targets.size(); k;) {
+                k -= 2;
+                buf->push_back(inst.targets[k]);
+                buf->push_back(inst.targets[k + 1]);
+            }
+            do_instruction_callback({GateType::CX, {}, *buf});
+        } else {
+            assert(inst.gate_type == GATE_DATA[inst.gate_type].inverse().id);
+            do_instruction_callback(inst);
+        }
+    });
+}
+
+void stim::decompose_cpp_operation_with_reverse_independence(
+    const CircuitInstruction &cpp_op,
+    size_t num_qubits,
+    const std::function<void(const CircuitInstruction &inst)> &do_instruction_callback) {
+    PauliString<64> obs1(num_qubits);
+    PauliString<64> obs2(num_qubits);
+    std::vector<GateTarget> bits1;
+    std::vector<GateTarget> bits2;
+    Circuit circuit_workspace;
+    std::vector<GateTarget> target_buf;
+
+    size_t start = 0;
+    while (true) {
+        bool b1 = accumulate_next_obs_terms_to_pauli_string_helper(cpp_op, &start, &obs1, &bits1);
+        bool b2 = accumulate_next_obs_terms_to_pauli_string_helper(cpp_op, &start, &obs2, &bits2);
+        if (!b2) {
+            break;
+        }
+        if (!b1) {
+            throw std::invalid_argument("Odd number of products.");
+        }
+
+        decompose_cpp_operation_with_reverse_independence_helper(
+            cpp_op, obs1, obs2, bits1, bits2, do_instruction_callback, &circuit_workspace, &target_buf);
+    }
+}
+
 void stim::decompose_pair_instruction_into_segments_with_single_use_controls(
     const CircuitInstruction &inst, size_t num_qubits, const std::function<void(CircuitInstruction)> &callback) {
     simd_bits<64> used_as_control(std::max(num_qubits, size_t{1}));
@@ -709,6 +859,11 @@ struct Simplifier {
                     simplify_instruction(sub);
                 });
                 break;
+            case GateType::CPP:
+                decompose_cpp_operation_with_reverse_independence(inst, num_qubits, [&](const CircuitInstruction sub) {
+                    simplify_instruction(sub);
+                });
+                break;
             case GateType::SPP:
             case GateType::SPP_DAG:
                 decompose_spp_or_spp_dag_operation(inst, num_qubits, false, [&](const CircuitInstruction sub) {

src/stim/circuit/gate_decomposition.h

@@ -59,6 +59,15 @@ void decompose_mpp_operation(
     size_t num_qubits,
     const std::function<void(const CircuitInstruction &inst)> &do_instruction_callback);
 
+/// Decomposes CPP operations into sequences of simpler operations with the same effect.
+///
+/// The output is guaranteed to only use self-inverse operations, and to have the same
+/// effect if run in order or in reverse order.
+void decompose_cpp_operation_with_reverse_independence(
+    const CircuitInstruction &cpp_op,
+    size_t num_qubits,
+    const std::function<void(const CircuitInstruction &inst)> &do_instruction_callback);
+
 /// Decomposes SPP operations into sequences of simpler operations with the same effect.
 void decompose_spp_or_spp_dag_operation(
     const CircuitInstruction &spp_op,