diff --git a/src/stim/stabilizers/conversions.h b/src/stim/stabilizers/conversions.h index 4349693fd..37166d70d 100644 --- a/src/stim/stabilizers/conversions.h +++ b/src/stim/stabilizers/conversions.h @@ -83,11 +83,13 @@ Circuit stabilizer_state_vector_to_circuit( /// ignore_noise: If the circuit contains noise channels, ignore them instead of raising an exception. /// ignore_measurement: If the circuit contains measurements, ignore them instead of raising an exception. /// ignore_reset: If the circuit contains resets, ignore them instead of raising an exception. +/// inverse: The last step of the implementation is to invert the tableau. Setting this argument +/// to true will skip this inversion, saving time but returning the inverse tableau. /// /// Returns: /// A tableau encoding the given circuit's Clifford operation. template -Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset); +Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset, bool inverse = false); /// Simulates the given circuit and outputs a state vector. /// diff --git a/src/stim/stabilizers/conversions.inl b/src/stim/stabilizers/conversions.inl index b6a8edec3..6ab5f2b17 100644 --- a/src/stim/stabilizers/conversions.inl +++ b/src/stim/stabilizers/conversions.inl @@ -149,7 +149,7 @@ std::vector>> tableau_to_unitary(const Tableau -Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset) { +Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset, bool inverse) { Tableau result(circuit.count_qubits()); TableauSimulator sim(std::mt19937_64(0), circuit.count_qubits()); @@ -185,7 +185,10 @@ Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ig } }); - return sim.inv_state.inverse(); + if (!inverse) { + return sim.inv_state.inverse(); + } + return sim.inv_state; } template @@ -556,7 +559,7 @@ Tableau stabilizers_to_tableau( } for (size_t k1 = 0; k1 < stabilizers.size(); k1++) { - for (size_t k2 = 0; k2 < stabilizers.size(); k2++) { + for (size_t k2 = k1 + 1; k2 < stabilizers.size(); k2++) { if (!stabilizers[k1].ref().commutes(stabilizers[k2])) { std::stringstream ss; ss << "Some of the given stabilizers anticommute.\n"; @@ -568,44 +571,39 @@ Tableau stabilizers_to_tableau( } } } - Tableau inverted(num_qubits); - - PauliString cur(num_qubits); - std::vector targets; - while (targets.size() < num_qubits) { - targets.push_back(targets.size()); - } - auto overwrite_cur_apply_recorded = [&](const PauliString &e) { - PauliStringRef cur_ref = cur.ref(); - cur.xs.clear(); - cur.zs.clear(); - cur.xs.word_range_ref(0, e.xs.num_simd_words) = e.xs; - cur.zs.word_range_ref(0, e.xs.num_simd_words) = e.zs; - cur.sign = e.sign; - inverted.apply_within(cur_ref, targets); - }; + Circuit elimination_instructions; + PauliString buf(num_qubits); size_t used = 0; for (const auto &e : stabilizers) { - overwrite_cur_apply_recorded(e); + if (e.num_qubits == num_qubits) { + buf = e; + } else { + buf.xs.clear(); + buf.zs.clear(); + memcpy(buf.xs.u8, e.xs.u8, e.xs.num_u8_padded()); + memcpy(buf.zs.u8, e.zs.u8, e.zs.num_u8_padded()); + buf.sign = e.sign; + } + buf.ref().do_circuit(elimination_instructions); // Find a non-identity term in the Pauli string past the region used by other stabilizers. size_t pivot; for (pivot = used; pivot < num_qubits; pivot++) { - if (cur.xs[pivot] || cur.zs[pivot]) { + if (buf.xs[pivot] || buf.zs[pivot]) { break; } } // Check for incompatible / redundant stabilizers. if (pivot == num_qubits) { - if (cur.xs.not_zero()) { + if (buf.xs.not_zero()) { throw std::invalid_argument("Some of the given stabilizers anticommute."); } - if (cur.sign) { + if (buf.sign) { throw std::invalid_argument("Some of the given stabilizers contradict each other."); } - if (!allow_redundant && cur.zs.not_zero()) { + if (!allow_redundant && buf.zs.not_zero()) { throw std::invalid_argument( "Didn't specify allow_redundant=True but one of the given stabilizers is a product of the others. " "To allow redundant stabilizers, pass the argument allow_redundant=True."); @@ -614,32 +612,36 @@ Tableau stabilizers_to_tableau( } // Change pivot basis to the Z axis. - if (cur.xs[pivot]) { - std::string name = cur.zs[pivot] ? "H_YZ" : "H_XZ"; - inverted.inplace_scatter_append(GATE_DATA.at(name).tableau(), {pivot}); + if (buf.xs[pivot]) { + GateType g = buf.zs[pivot] ? GateType::H_YZ : GateType::H; + GateTarget t = GateTarget::qubit(pivot); + CircuitInstruction instruction{g, {}, &t}; + elimination_instructions.safe_append(instruction); + buf.ref().do_instruction(instruction); } // Cancel other terms in Pauli string. for (size_t q = 0; q < num_qubits; q++) { - int p = cur.xs[q] + cur.zs[q] * 2; + int p = buf.xs[q] + buf.zs[q] * 2; if (p && q != pivot) { - inverted.inplace_scatter_append( - GATE_DATA.at(p == 1 ? "XCX" - : p == 2 ? "XCZ" - : "XCY") - .tableau(), - {pivot, q}); + std::array targets{GateTarget::qubit(pivot), GateTarget::qubit(q)}; + CircuitInstruction instruction{p == 1 ? GateType::XCX : p == 2 ? GateType::XCZ : GateType::XCY, {}, targets}; + elimination_instructions.safe_append(instruction); + buf.ref().do_instruction(instruction); } } // Move pivot to diagonal. if (pivot != used) { - inverted.inplace_scatter_append(GATE_DATA.at("SWAP").tableau(), {pivot, used}); + std::array targets{GateTarget::qubit(pivot), GateTarget::qubit(used)}; + CircuitInstruction instruction{GateType::SWAP, {}, targets}; + elimination_instructions.safe_append(instruction); } // Fix sign. - overwrite_cur_apply_recorded(e); - if (cur.sign) { - inverted.inplace_scatter_append(GATE_DATA.at("X").tableau(), {used}); + if (buf.sign) { + GateTarget t = GateTarget::qubit(used); + CircuitInstruction instruction{GateType::X, {}, &t}; + elimination_instructions.safe_append(instruction); } used++; @@ -653,10 +655,17 @@ Tableau stabilizers_to_tableau( } } + if (num_qubits > 0) { + // Force size of resulting tableau to be correct. + GateTarget t = GateTarget::qubit(num_qubits - 1); + elimination_instructions.safe_append(CircuitInstruction{GateType::X, {}, &t}); + elimination_instructions.safe_append(CircuitInstruction{GateType::X, {}, &t}); + } + if (invert) { - return inverted; + return circuit_to_tableau(elimination_instructions.inverse(), false, false, false, true); } - return inverted.inverse(); + return circuit_to_tableau(elimination_instructions, false, false, false, true); } } // namespace stim diff --git a/src/stim/stabilizers/conversions.perf.cc b/src/stim/stabilizers/conversions.perf.cc index 0461355a7..50ad12c50 100644 --- a/src/stim/stabilizers/conversions.perf.cc +++ b/src/stim/stabilizers/conversions.perf.cc @@ -67,3 +67,51 @@ BENCHMARK(independent_to_disjoint_xyz_errors) { std::cout << "data dependence"; } } + +BENCHMARK(stabilizers_to_tableau) { + std::vector> offsets{ + {1, 0}, + {-1, 0}, + {0, 1}, + {0, -1}, + {3, 6}, + {-6, 3}, + }; + size_t w = 24; + size_t h = 12; + + auto normalize = [&](std::complex c) -> std::complex { + return {fmodf(c.real() + w*10, w), fmodf(c.imag() + h*10, h)}; + }; + auto q2i = [&](std::complex c) -> size_t { + c = normalize(c); + return (int)c.real() / 2 + c.imag() * (w / 2); + }; + + std::vector> stabilizers; + for (size_t x = 0; x < w; x++) { + for (size_t y = x % 2; y < h; y += 2) { + std::complex s{x % 2 ? -1.0f : +1.0f, 0.0f}; + std::complex c{(float)x, (float)y}; + stim::PauliString<64> ps(w * h / 2); + for (const auto &offset : offsets) { + size_t i = q2i(c + offset * s); + if (x % 2 == 0) { + ps.xs[i] = 1; + } else { + ps.zs[i] = 1; + } + } + stabilizers.push_back(ps); + } + } + + size_t dep = 0; + benchmark_go([&]() { + Tableau<64> t = stabilizers_to_tableau(stabilizers, true, true, false); + dep += t.xs[0].zs[0]; + }).goal_millis(5); + if (dep == 99999999) { + std::cout << "data dependence"; + } +}