diff --git a/Cargo.toml b/Cargo.toml index 805450ac..51dc243b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,6 +5,7 @@ members = [ "paulimer", "paulimer/bindings/python", "pauliverse", + "test-utils/dense-oracle", "deq/deq_runtime", "deq/deq_decoder_abi", "deq/deq_decoder_abi/reference_plugin", diff --git a/README.md b/README.md index af76f05b..d75a91a4 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ This repository contains several interconnected crates: - [binar](binar): A high-performance bit manipulation library providing bit vectors, bit matrices, and bitwise operations over GF(2). - [paulimer](paulimer): A library for Pauli operators and Clifford gates, built on binar. -- [pauliverse](pauliverse): Fast stabilizer simulators. +- [pauliverse](pauliverse): Fast stabilizer simulators, including exact global-phase tracking (`PhasedOutcomeCompleteSimulation`) for verifying parameterised circuits. - [deq](deq): A dynamic and generic QEC decoding system, including the `.deq` DSL, transpiler, JIT runtime (Rust), CLI, and an anywidget-based visualizer. ### Python Bindings @@ -136,6 +136,26 @@ let image = clifford.image(&x0); assert_eq!(image, "XX".parse::().unwrap()); ``` +### Stabilizer Simulation (pauliverse) + +```rust +use pauliverse::{PhasedOutcomeCompleteSimulation, Simulation}; +use paulimer::{SparsePauli, UnitaryOp}; + +// Track the exact global phase while simulating all measurement-outcome branches. +let mut sim = PhasedOutcomeCompleteSimulation::new(2); +sim.unitary_op(UnitaryOp::Hadamard, &[0]); +sim.unitary_op(UnitaryOp::Hadamard, &[1]); + +// Apply a symbolic rotation e^{i·alpha·Z0Z1} around a higher-weight (two-qubit) Pauli. +let alpha = sim.allocate_symbolic_angle(); +sim.symbolic_pauli_exp(&"Z0 Z1".parse::().unwrap(), alpha); +``` + +The phase-aware `phased_action` of two circuits can then be compared for exact equality — see the +[pauliverse README](pauliverse/README.md) and the Python +[example notebooks](paulimer/bindings/python/examples). + ## Benchmarks This repository uses [Criterion](https://github.com/bheisler/criterion.rs) for Rust benchmarks and [ASV](https://asv.readthedocs.io/) for Python benchmarks. diff --git a/paulimer/Cargo.toml b/paulimer/Cargo.toml index 7d6fc2b8..c429c697 100644 --- a/paulimer/Cargo.toml +++ b/paulimer/Cargo.toml @@ -39,6 +39,7 @@ criterion = { version = "0.5", features = ["html_reports"] } proptest = "1.0" serde_json = "1.0" jsonschema = { version = "0.29" } +dense-oracle = { path = "../test-utils/dense-oracle" } [[bench]] name = "pauli_benchmark" diff --git a/paulimer/README.md b/paulimer/README.md index 66215631..e72393a1 100644 --- a/paulimer/README.md +++ b/paulimer/README.md @@ -21,6 +21,9 @@ the building blocks for stabilizer quantum mechanics and quantum error correctio - **Clifford Unitaries**: Efficient representation enabling fast operations - [`CliffordUnitary`]: O(n²) Pauli conjugation via binary symplectic matrix - Supports all standard Clifford gates (H, S, CNOT, etc.) + - [`clifford_to_pauli_exponents`]: exact decomposition into `π/4` Pauli exponents (the full + tableau, including image signs), so replaying it with exact phase tracking yields a well-defined + global phase Based on algorithms from [arXiv:2309.08676](https://arxiv.org/abs/2309.08676). @@ -94,6 +97,15 @@ assert_eq!(image, "XX".parse::().unwrap()); let mut circuit = CliffordUnitary::identity(2); circuit.left_mul(UnitaryOp::Hadamard, &[0]); circuit.left_mul(UnitaryOp::ControlledX, &[0, 1]); + +// Decompose a Clifford into an ordered product of π/4 Pauli exponents (exact, sign-preserving) +use paulimer::clifford::clifford_to_pauli_exponents; +let exponents = clifford_to_pauli_exponents(&circuit); +let mut rebuilt = CliffordUnitary::identity(2); +for pauli in &exponents { + rebuilt.left_mul_pauli_exp(pauli); +} +assert_eq!(rebuilt, circuit); ``` ## When to Use Each Type @@ -171,6 +183,7 @@ Key documentation: - [`SparsePauli`](src/pauli/sparse.rs) - Sparse Pauli representation for large systems - [`PauliGroup`](src/pauli_group.rs) - Subgroup operations and stabilizer groups - [`CliffordUnitary`](src/clifford.rs) - Clifford gates and Pauli conjugation +- [`clifford_to_pauli_exponents`](src/clifford/decomposition.rs) - Exact decomposition of a Clifford into `π/4` Pauli exponents - [Trait documentation](src/lib.rs) - `Pauli`, `Clifford`, and other core traits ## Contributing diff --git a/paulimer/bindings/python/README.md b/paulimer/bindings/python/README.md index 281d240c..3f9e0199 100644 --- a/paulimer/bindings/python/README.md +++ b/paulimer/bindings/python/README.md @@ -29,16 +29,71 @@ sim.apply_unitary(paulimer.UnitaryOpcode.ControlledX, [0, 1]) sim.measure(paulimer.SparsePauli("Z0")) ``` +### Verifying parameterised circuits with symbolic angles + +`PhasedOutcomeCompleteSimulation` tracks the exact global phase, so two circuits that share the same +free rotation angles can be checked for exact equality — even for exponents of higher-weight Paulis: + +```python +from paulimer import PhasedOutcomeCompleteSimulation, SparsePauli, UnitaryOpcode + + +def prepared_action(build): + sim = PhasedOutcomeCompleteSimulation(2) + for qubit in range(2): + sim.apply_unitary(UnitaryOpcode.Hadamard, [qubit]) # |++> + build(sim) + return sim.phased_action([], [0, 1]) # state-preparation action + + +def direct(sim): # e^{i alpha Z0 Z1} |++> + alpha = sim.allocate_symbolic_angle() + sim.apply_symbolic_pauli_exp(SparsePauli("Z_0 Z_1"), alpha) + + +def conjugated(sim): # CNOT . e^{i alpha Z1} . CNOT |++> + sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + alpha = sim.allocate_symbolic_angle() + sim.apply_symbolic_pauli_exp(SparsePauli("Z_1"), alpha) + sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + + +assert prepared_action(direct).is_equivalent(prepared_action(conjugated)) +``` + +### Decomposing a Clifford into pi/4 Pauli exponents + +`CliffordUnitary.to_pauli_exponents()` returns an ordered product of `pi/4` Pauli exponents that +reproduces the Clifford exactly, including the Pauli-image signs — so replaying it with exact phase +tracking yields a well-defined global phase: + +```python +from paulimer import CliffordUnitary, UnitaryOpcode + +clifford = CliffordUnitary.identity(2) +clifford.left_mul(UnitaryOpcode.Hadamard, [0]) +clifford.left_mul(UnitaryOpcode.ControlledX, [0, 1]) + +exponents = clifford.to_pauli_exponents() # list[SparsePauli], each factor exp(+-i pi/4 P) + +rebuilt = CliffordUnitary.identity(2) +for pauli in exponents: + rebuilt.left_mul_pauli_exp(pauli) +assert rebuilt == clifford +``` + ## Features - **DensePauli / SparsePauli** - Pauli operators with phase tracking and multiplication - **CliffordUnitary** - Clifford gates with conjugation and composition + (including `to_pauli_exponents()`, an exact decomposition into `pi/4` Pauli exponents) - **PauliGroup** - Group operations including membership testing and factorization - **Stabilizer Simulation** - Noiseless (OutcomeComplete, OutcomeFree, OutcomeSpecific) and noisy (Faulty) modes +- **PhasedOutcomeCompleteSimulation** - Outcome-complete simulation that additionally tracks the exact global phase, enabling exact equality checking of parameterised (symbolic-angle) circuits ## Use Cases -Designed for quantum error correction research, including stabilizer circuit analysis and Clifford circuit verification. +Designed for quantum error correction research, including stabilizer circuit analysis and Clifford circuit verification. `PhasedOutcomeCompleteSimulation` extends this to exact verification of non-stabilizer circuits built from symbolic Pauli rotations `e^{i alpha P}`. ## Performance diff --git a/paulimer/bindings/python/examples/phased-outcome-complete-simulation.ipynb b/paulimer/bindings/python/examples/phased-outcome-complete-simulation.ipynb new file mode 100644 index 00000000..f92acb03 --- /dev/null +++ b/paulimer/bindings/python/examples/phased-outcome-complete-simulation.ipynb @@ -0,0 +1,342 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6ab8d189", + "metadata": {}, + "source": [ + "# Tracking the Exact Global Phase with `PhasedOutcomeCompleteSimulation`\n", + "\n", + "This notebook introduces `PhasedOutcomeCompleteSimulation`, the global-phase-resolving\n", + "generalization of `OutcomeCompleteSimulation`. It implements Algorithm 4.2 of\n", + "[arXiv:2603.24717](https://arxiv.org/abs/2603.24717), which extends the outcome-complete\n", + "stabilizer simulation of [arXiv:2309.08676](https://arxiv.org/abs/2309.08676) to keep track of\n", + "the **exact global phase** of the simulated state.\n", + "\n", + "## Why the global phase matters\n", + "\n", + "Ordinary stabilizer simulation only represents a state *up to a global phase*. That is enough for\n", + "sampling measurement outcomes, but **not** enough to verify equivalence of non-stabilizer circuits.\n", + "The motivating example from the paper is checking, for two pairs of Clifford unitaries\n", + "$(C_1, C_2)$ and $(D_1, D_2)$, whether\n", + "\n", + "$$ C_1\\, e^{i\\alpha Z_1}\\, C_2\\,|0\\rangle \\;=\\; D_1\\, e^{i\\alpha Z_1}\\, D_2\\,|0\\rangle \\quad\\text{for all } \\alpha. $$\n", + "\n", + "This holds **iff** the stabilizer states $C_1 Z_1^a C_2 |0\\rangle$ and $D_1 Z_1^a D_2 |0\\rangle$ are\n", + "equal for $a \\in \\{0, 1\\}$ — and that equality must be **exact, including the global phase**.\n", + "`PhasedOutcomeCompleteSimulation` is the engine that makes this exact comparison possible.\n", + "\n", + "Crucially, none of this ever forms an exponentially large state vector: phases are tracked\n", + "symbolically as integer powers of $\\zeta_8 = e^{i\\pi/4}$, and every operation stays polynomial in the\n", + "number of qubits." + ] + }, + { + "cell_type": "markdown", + "id": "657772da", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "baf07efc", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-27T19:52:07.173301Z", + "iopub.status.busy": "2026-06-27T19:52:07.173182Z", + "iopub.status.idle": "2026-06-27T19:52:07.176394Z", + "shell.execute_reply": "2026-06-27T19:52:07.175818Z" + } + }, + "outputs": [], + "source": [ + "from paulimer import (\n", + " PhasedOutcomeCompleteSimulation,\n", + " OutcomeCompleteSimulation,\n", + " SparsePauli,\n", + " UnitaryOpcode,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "da6db538", + "metadata": {}, + "source": [ + "## The phased stabilizer-state representation\n", + "\n", + "Like `OutcomeCompleteSimulation`, the phased simulator tracks **all** $2^{n_r}$ branches of a\n", + "circuit with $n_r$ random measurement outcomes at once. For a random-bit assignment\n", + "$r \\in \\{0,1\\}^{n_r}$ the encoded state is\n", + "\n", + "$$ i^{\\langle p,\\, r\\rangle}\\, (-1)^{\\langle B r + s,\\, r\\rangle}\\; R\\,|A r\\rangle, $$\n", + "\n", + "where\n", + "\n", + "- `R` is the phased state encoder (a Clifford unitary *with* its exact global phase),\n", + "- `A` = `sign_matrix` maps the random bits to a computational-basis input of `R`,\n", + "- `B` = `quadratic_phase_matrix` carries the quadratic $(-1)$-phase,\n", + "- `p` = `linear_i_phase` carries the linear $i$-phase,\n", + "- `s` = `linear_sign_phase` carries the linear $(-1)$-phase.\n", + "\n", + "The scalar prefactor $i^{\\langle p,r\\rangle}(-1)^{\\langle Br+s,r\\rangle}$ is returned, as a\n", + "$\\zeta_8$ exponent in $\\{0,\\dots,7\\}$, by `output_phase_exponent(r)`." + ] + }, + { + "cell_type": "markdown", + "id": "d8af2563", + "metadata": {}, + "source": [ + "## A worked example: interfering measurement phases\n", + "\n", + "We prepare $|{+}{+}\\rangle$ and measure the Pauli $Y$ on each qubit. Each $Y$ measurement is a random\n", + "outcome, so there are two random bits $r = (r_0, r_1)$ and four branches." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a4a409ce", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-27T19:52:07.178916Z", + "iopub.status.busy": "2026-06-27T19:52:07.178852Z", + "iopub.status.idle": "2026-06-27T19:52:07.181425Z", + "shell.execute_reply": "2026-06-27T19:52:07.181152Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "qubits : 2\n", + "random outcomes : 2\n", + "A (sign_matrix) : 10\n", + "01\n", + "\n", + "B (quad. phase) : 01\n", + "00\n", + "\n", + "p (linear i) : [11]\n", + "s (linear sign) : [00]\n" + ] + } + ], + "source": [ + "sim = PhasedOutcomeCompleteSimulation(2)\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0])\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [1])\n", + "sim.measure(SparsePauli(\"Y_0\"))\n", + "sim.measure(SparsePauli(\"Y_1\"))\n", + "\n", + "print(\"qubits :\", sim.qubit_count)\n", + "print(\"random outcomes :\", sim.random_outcome_count)\n", + "print(\"A (sign_matrix) :\", repr(sim.sign_matrix))\n", + "print(\"B (quad. phase) :\", repr(sim.quadratic_phase_matrix))\n", + "print(\"p (linear i) :\", sim.linear_i_phase)\n", + "print(\"s (linear sign) :\", sim.linear_sign_phase)" + ] + }, + { + "cell_type": "markdown", + "id": "be6d0508", + "metadata": {}, + "source": [ + "The linear $i$-phase is $p = (1, 1)$: each random $Y$ outcome contributes a factor of $i$.\n", + "The quadratic matrix $B$ has a single off-diagonal entry. Let us read off the $\\zeta_8$ exponent of\n", + "every branch and translate it into a familiar scalar." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "67772ee4", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-27T19:52:07.182817Z", + "iopub.status.busy": "2026-06-27T19:52:07.182729Z", + "iopub.status.idle": "2026-06-27T19:52:07.184643Z", + "shell.execute_reply": "2026-06-27T19:52:07.184505Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "r = (0, 0) -> zeta8^0 = 1\n", + "r = (1, 0) -> zeta8^2 = i\n", + "r = (0, 1) -> zeta8^2 = i\n", + "r = (1, 1) -> zeta8^4 = -1\n" + ] + } + ], + "source": [ + "ZETA8_LABEL = {0: \"1\", 1: \"ζ₈\", 2: \"i\", 3: \"ζ₈³\", 4: \"-1\", 5: \"ζ₈⁵\", 6: \"-i\", 7: \"ζ₈⁷\"}\n", + "\n", + "def branches(n_random):\n", + " for k in range(2 ** n_random):\n", + " yield [bool((k >> i) & 1) for i in range(n_random)]\n", + "\n", + "for r in branches(sim.random_outcome_count):\n", + " e = sim.output_phase_exponent(r)\n", + " print(f\"r = {tuple(int(b) for b in r)} -> zeta8^{e} = {ZETA8_LABEL[e]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "048862d9", + "metadata": {}, + "source": [ + "## The punchline\n", + "\n", + "Look at the four phases: $1,\\ i,\\ i,\\ -1$.\n", + "\n", + "- Branches $(1,0)$ and $(0,1)$ each pick up a single $i$ from one $Y$ outcome.\n", + "- Branch $(1,1)$ picks up **both**, and the two factors of $i$ interfere to give\n", + " $i \\cdot i = -1$ (that is $\\zeta_8^4$), *not* $i + i$.\n", + "\n", + "This is exactly the kind of exact, non-trivial global phase that ordinary stabilizer simulation\n", + "discards. The phased simulator stores the linear $i$-phase $p$ **modulo 2**, so the $i\\cdot i = -1$\n", + "carry is absorbed into the quadratic matrix $B$ — which is why $B$ is needed at all.\n", + "\n", + "The plain `OutcomeCompleteSimulation` simply does not expose this information:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "6e37913e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-27T19:52:07.187224Z", + "iopub.status.busy": "2026-06-27T19:52:07.187174Z", + "iopub.status.idle": "2026-06-27T19:52:07.188648Z", + "shell.execute_reply": "2026-06-27T19:52:07.188376Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "OutcomeCompleteSimulation has output_phase_exponent? False\n", + "OutcomeCompleteSimulation has linear_i_phase? False\n" + ] + } + ], + "source": [ + "plain = OutcomeCompleteSimulation(2)\n", + "print(\"OutcomeCompleteSimulation has output_phase_exponent? \",\n", + " hasattr(plain, \"output_phase_exponent\"))\n", + "print(\"OutcomeCompleteSimulation has linear_i_phase? \",\n", + " hasattr(plain, \"linear_i_phase\"))" + ] + }, + { + "cell_type": "markdown", + "id": "1ffdb295", + "metadata": {}, + "source": [ + "## Recomputing the phase from the raw data\n", + "\n", + "To make the formula concrete, we recompute $i^{\\langle p,r\\rangle}(-1)^{\\langle Br+s,r\\rangle}$\n", + "directly from the exposed matrices and vectors — using only integer (GF(2)) arithmetic — and check\n", + "it against `output_phase_exponent` for every branch." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "031ad2f8", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-27T19:52:07.189602Z", + "iopub.status.busy": "2026-06-27T19:52:07.189555Z", + "iopub.status.idle": "2026-06-27T19:52:07.191532Z", + "shell.execute_reply": "2026-06-27T19:52:07.191389Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "✓ Reconstructed phases match output_phase_exponent for all branches\n" + ] + } + ], + "source": [ + "def phase_exponent_from_data(sim, r):\n", + " \"\"\"Reconstruct the zeta8 exponent from A, B, p, s using exact integer arithmetic.\"\"\"\n", + " B = sim.quadratic_phase_matrix\n", + " p = sim.linear_i_phase\n", + " s = sim.linear_sign_phase\n", + " n = sim.random_outcome_count\n", + "\n", + " linear_i = sum(1 for i in range(n) if p[i] and r[i]) % 2\n", + "\n", + " quadratic = 0\n", + " for i in range(n):\n", + " for j in range(n):\n", + " if B[i, j] and r[i] and r[j]:\n", + " quadratic ^= 1\n", + " linear_sign = sum(1 for i in range(n) if s[i] and r[i]) % 2\n", + " sign = quadratic ^ linear_sign\n", + "\n", + " return (2 * linear_i + 4 * sign) % 8\n", + "\n", + "for r in branches(sim.random_outcome_count):\n", + " assert phase_exponent_from_data(sim, r) == sim.output_phase_exponent(r)\n", + "print(\"✓ Reconstructed phases match output_phase_exponent for all branches\")" + ] + }, + { + "cell_type": "markdown", + "id": "f2fdfbbe", + "metadata": {}, + "source": [ + "## Why this enables circuit verification\n", + "\n", + "Recall the equivalence $C_1 e^{i\\alpha Z_1} C_2 |0\\rangle = D_1 e^{i\\alpha Z_1} D_2|0\\rangle$\n", + "reduces to the **exact** equality of the stabilizer states $C_1 Z_1^a C_2 |0\\rangle$ and\n", + "$D_1 Z_1^a D_2 |0\\rangle$ for $a \\in \\{0,1\\}$. Because `PhasedOutcomeCompleteSimulation` retains the\n", + "exact global phase of each branch's encoded state $R|Ar\\rangle$ — rather than only the state up to a\n", + "phase — it provides precisely the information needed to decide such equalities. The same machinery\n", + "extends to circuits with intermediate measurements and outcome-parity-conditional Pauli gates, which\n", + "are ubiquitous in fault-tolerant quantum computing.\n", + "\n", + "Everything above runs in time polynomial in the number of qubits: the simulator never materializes a\n", + "$2^n$-dimensional state vector, and the global phase is carried exactly as an integer $\\zeta_8$\n", + "exponent. The implementation is validated phase-exactly against brute-force dense statevector\n", + "references in the Rust test suite (`pauliverse/tests/phased_outcome_complete_dense.rs`)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/paulimer/bindings/python/examples/verifying-circuit-equivalence.ipynb b/paulimer/bindings/python/examples/verifying-circuit-equivalence.ipynb new file mode 100644 index 00000000..d50bb561 --- /dev/null +++ b/paulimer/bindings/python/examples/verifying-circuit-equivalence.ipynb @@ -0,0 +1,304 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "53409630", + "metadata": {}, + "source": [ + "# Verifying Circuit Equivalence (Section 4.1)\n", + "\n", + "This notebook illustrates the **circuit-equivalence verification** application of\n", + "[arXiv:2603.24717](https://arxiv.org/abs/2603.24717), Section 4.1, using the Python\n", + "`PhasedOutcomeCompleteSimulation` API. The companion notebook *Verifying Symbolic Pauli-Exponent Circuits*\n", + "covers the same idea for **channels** (via Choi states); here we focus on the literal Section 4.1\n", + "construction for **state preparation**.\n", + "\n", + "## The reduction\n", + "\n", + "We want to decide whether two parameterized state-preparation circuits agree for *every* angle $\\alpha$:\n", + "$$ C_1\\, e^{i\\alpha Z}\\, C_2\\,|0\\cdots0\\rangle \\;=\\; D_1\\, e^{i\\alpha Z}\\, D_2\\,|0\\cdots0\\rangle \\quad\\text{for all }\\alpha. $$\n", + "Section 4.1 shows this is equivalent to a *single* **exact** stabilizer-state equality, with the\n", + "continuous angle replaced by a binary symbolic exponent $a$:\n", + "$$ C_1\\, Z^{a}\\, C_2\\,|0\\cdots0\\rangle \\;=\\; D_1\\, Z^{a}\\, D_2\\,|0\\cdots0\\rangle. $$\n", + "Exactness is the whole point: the equality must hold *including* the relative phase between the\n", + "$a=0$ and $a=1$ branches, since that phase is what encodes the Pauli exponent for every $\\alpha$. An\n", + "ordinary, phase-blind stabilizer comparison is not enough.\n", + "\n", + "There is **no dedicated verification function**: the check is simply `phased_action(...)` followed by\n", + "`PhasedCircuitAction.is_equivalent(...)`, exactly mirroring how `OutcomeCompleteSimulation` performs\n", + "phaseless equality checking. The symbolic exponent $Z^{a}$ is added with `allocate_symbolic_angle()` and\n", + "`apply_symbolic_pauli_exp`." + ] + }, + { + "cell_type": "markdown", + "id": "3286f077", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "63502173", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-28T03:55:25.320063Z", + "iopub.status.busy": "2026-06-28T03:55:25.319936Z", + "iopub.status.idle": "2026-06-28T03:55:25.325161Z", + "shell.execute_reply": "2026-06-28T03:55:25.324077Z" + } + }, + "outputs": [], + "source": [ + "from paulimer import (\n", + " PhasedOutcomeCompleteSimulation,\n", + " PhasedCircuitAction,\n", + " SparsePauli,\n", + " UnitaryOpcode,\n", + ")\n", + "\n", + "\n", + "def prepared_state_action(qubit_count, build):\n", + " \"\"\"Record a state-preparation circuit C_1 (prod_k Z_k^{a_k}) C_2 |0...0> as a phased action.\n", + "\n", + " State preparation is the `inputs = []` case of `phased_action`; the outputs are all qubits.\n", + " \"\"\"\n", + " sim = PhasedOutcomeCompleteSimulation(qubit_count)\n", + " build(sim)\n", + " return sim.phased_action([], list(range(qubit_count)))\n", + "\n", + "\n", + "def prepare_plus(sim, qubit_count):\n", + " \"\"\"Prepare |+...+> by Hadamarding every qubit.\"\"\"\n", + " for qubit in range(qubit_count):\n", + " sim.apply_unitary(UnitaryOpcode.Hadamard, [qubit])" + ] + }, + { + "cell_type": "markdown", + "id": "f85ba10d", + "metadata": {}, + "source": [ + "## Two factorizations of the same state\n", + "\n", + "A clean Section 4.1 instance: the parameterized state $e^{i\\alpha Z_0 Z_1}\\,|{+}{+}\\rangle$ written\n", + "two different ways. Conjugating a single-qubit Pauli exponent by a `CNOT` turns it into a two-qubit $ZZ$\n", + "exponent, because $\\mathrm{CNOT}_{01}\\,Z_1\\,\\mathrm{CNOT}_{01} = Z_0 Z_1$ and\n", + "$\\mathrm{CNOT}_{01}\\,|{+}{+}\\rangle = |{+}{+}\\rangle$:\n", + "$$ e^{i\\alpha Z_0 Z_1}\\,|{+}{+}\\rangle \\;=\\; \\mathrm{CNOT}_{01}\\; e^{i\\alpha Z_1}\\; \\mathrm{CNOT}_{01}\\,|{+}{+}\\rangle. $$\n", + "The two circuits are different Clifford factorizations $C_1 Z^a C_2$ of the same parameterized state,\n", + "so they must verify as equivalent." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "864afa8f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-28T03:55:25.327371Z", + "iopub.status.busy": "2026-06-28T03:55:25.327314Z", + "iopub.status.idle": "2026-06-28T03:55:25.330312Z", + "shell.execute_reply": "2026-06-28T03:55:25.329545Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified equal: e^{i alpha Z0Z1}|++> == CNOT01 . e^{i alpha Z1} . CNOT01 |++>\n" + ] + } + ], + "source": [ + "def direct_zz(sim): # e^{i alpha Z0Z1} |++> (C_2 = H(x)H, Z^a on Z0Z1, C_1 = I)\n", + " prepare_plus(sim, 2)\n", + " alpha = sim.allocate_symbolic_angle()\n", + " sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0 Z_1\"), alpha)\n", + "\n", + "\n", + "def conjugated_zz(sim): # CNOT01 . e^{i alpha Z1} . CNOT01 |++>\n", + " prepare_plus(sim, 2)\n", + " sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + " alpha = sim.allocate_symbolic_angle()\n", + " sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_1\"), alpha)\n", + " sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + "\n", + "\n", + "direct = prepared_state_action(2, direct_zz)\n", + "conjugated = prepared_state_action(2, conjugated_zz)\n", + "\n", + "assert direct.is_equivalent(conjugated)\n", + "assert conjugated.is_equivalent(direct) # verification is symmetric\n", + "print(\"verified equal: e^{i alpha Z0Z1}|++> == CNOT01 . e^{i alpha Z1} . CNOT01 |++>\")" + ] + }, + { + "cell_type": "markdown", + "id": "96811100", + "metadata": {}, + "source": [ + "## A purely-phase difference\n", + "\n", + "The verification is phase-sensitive. The states $e^{+i\\alpha Z}\\,|{+}\\rangle$ and\n", + "$e^{-i\\alpha Z}\\,|{+}\\rangle$ have **identical** stabilizer (symplectic) data, differing only by a\n", + "relative $-1$ on the $a=1$ branch. A phase-blind comparison (`is_equivalent_up_to_signs`) reports them\n", + "equal; the phase-aware `is_equivalent` correctly separates them." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6a34aaf7", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-28T03:55:25.332317Z", + "iopub.status.busy": "2026-06-28T03:55:25.332271Z", + "iopub.status.idle": "2026-06-28T03:55:25.336344Z", + "shell.execute_reply": "2026-06-28T03:55:25.334863Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "phase-blind: e^{+i alpha Z}|+> and e^{-i alpha Z}|+> are INDISTINGUISHABLE\n", + "phase-aware: e^{+i alpha Z}|+> != e^{-i alpha Z}|+>\n" + ] + } + ], + "source": [ + "def rotate_plus(sign):\n", + " pauli = SparsePauli(\"Z_0\") if sign > 0 else SparsePauli(\"-Z_0\")\n", + "\n", + " def build(sim):\n", + " prepare_plus(sim, 1)\n", + " alpha = sim.allocate_symbolic_angle()\n", + " sim.apply_symbolic_pauli_exp(pauli, alpha)\n", + "\n", + " return prepared_state_action(1, build)\n", + "\n", + "\n", + "positive, negative = rotate_plus(+1), rotate_plus(-1)\n", + "\n", + "assert positive.is_equivalent_up_to_signs(negative) # phase-blind: indistinguishable\n", + "assert not positive.is_equivalent(negative) # phase-aware: distinct\n", + "print(\"phase-blind: e^{+i alpha Z}|+> and e^{-i alpha Z}|+> are INDISTINGUISHABLE\")\n", + "print(\"phase-aware: e^{+i alpha Z}|+> != e^{-i alpha Z}|+>\")" + ] + }, + { + "cell_type": "markdown", + "id": "40f0f2f9", + "metadata": {}, + "source": [ + "## Several independent angles\n", + "\n", + "The reduction generalizes to several independent symbolic angles at once. We verify\n", + "$e^{i\\alpha Z_0 Z_1}\\, e^{i\\beta Z_0}\\,|{+}{+}\\rangle$ against its CNOT-conjugated factorization. The\n", + "two angles are allocated in the **same order** in both circuits, so the symbolic angles correspond\n", + "one-to-one. Negating the second exponent's Pauli injects a pure branch-phase difference, which is\n", + "detected." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "82542c1c", + "metadata": { + "execution": { + "iopub.execute_input": "2026-06-28T03:55:25.337755Z", + "iopub.status.busy": "2026-06-28T03:55:25.337671Z", + "iopub.status.idle": "2026-06-28T03:55:25.342463Z", + "shell.execute_reply": "2026-06-28T03:55:25.340361Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified equal: e^{i alpha Z0Z1} e^{i beta Z0}|++> == CNOT-conjugated factorization\n", + "detected: negating the second exponent breaks the branch-phase equality\n" + ] + } + ], + "source": [ + "def two_angle_direct(negate_second):\n", + " second = SparsePauli(\"-Z_0\") if negate_second else SparsePauli(\"Z_0\")\n", + "\n", + " def build(sim):\n", + " prepare_plus(sim, 2)\n", + " alpha = sim.allocate_symbolic_angle() # alpha (allocated first)\n", + " sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0 Z_1\"), alpha)\n", + " beta = sim.allocate_symbolic_angle() # beta (allocated second)\n", + " sim.apply_symbolic_pauli_exp(second, beta)\n", + "\n", + " return prepared_state_action(2, build)\n", + "\n", + "\n", + "def two_angle_conjugated(sim):\n", + " prepare_plus(sim, 2)\n", + " sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + " alpha = sim.allocate_symbolic_angle() # alpha\n", + " sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_1\"), alpha)\n", + " sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + " beta = sim.allocate_symbolic_angle() # beta\n", + " sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0\"), beta)\n", + "\n", + "\n", + "conjugated = prepared_state_action(2, two_angle_conjugated)\n", + "\n", + "assert two_angle_direct(negate_second=False).is_equivalent(conjugated)\n", + "print(\"verified equal: e^{i alpha Z0Z1} e^{i beta Z0}|++> == CNOT-conjugated factorization\")\n", + "\n", + "assert not two_angle_direct(negate_second=True).is_equivalent(conjugated)\n", + "print(\"detected: negating the second exponent breaks the branch-phase equality\")" + ] + }, + { + "cell_type": "markdown", + "id": "20d43148", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "- Section 4.1 reduces *all-angles* equivalence of $C_1 e^{i\\alpha Z} C_2|0\\rangle$ and\n", + " $D_1 e^{i\\alpha Z} D_2|0\\rangle$ to a single **exact** stabilizer-state equality\n", + " $C_1 Z^a C_2|0\\rangle = D_1 Z^a D_2|0\\rangle$, with $Z^a$ added via `allocate_symbolic_angle()` and `apply_symbolic_pauli_exp`.\n", + "\n", + "- The equality check needs **no special API**: build each state's `phased_action([], outputs)` and\n", + " compare with `PhasedCircuitAction.is_equivalent` (phase-aware) or `is_equivalent_up_to_signs`\n", + " (phase-blind). This is the phased counterpart of how `OutcomeCompleteSimulation` checks phaseless\n", + " equivalence.\n", + "- Because the comparison is **exact**, it certifies the relative branch phase, which is precisely what\n", + " distinguishes $e^{+i\\alpha Z}$ from $e^{-i\\alpha Z}$ and what ordinary stabilizer simulation discards." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/paulimer/bindings/python/examples/verifying-symbolic-rotations.ipynb b/paulimer/bindings/python/examples/verifying-symbolic-rotations.ipynb new file mode 100644 index 00000000..d679ef35 --- /dev/null +++ b/paulimer/bindings/python/examples/verifying-symbolic-rotations.ipynb @@ -0,0 +1,572 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f8907e83", + "metadata": {}, + "source": [ + "# Verifying Symbolic Pauli-Exponent Circuits with `PhasedOutcomeCompleteSimulation`\n", + "\n", + "This notebook shows how to apply a **symbolic Pauli exponent** $e^{i\\alpha P}$ to a\n", + "`PhasedOutcomeCompleteSimulation` and use it to verify equivalence of non-stabilizer circuits, the\n", + "motivating application of [arXiv:2603.24717](https://arxiv.org/abs/2603.24717). (See the companion\n", + "notebook *Tracking the Exact Global Phase* for an introduction to the phased simulator itself.)\n", + "\n", + "## Applying a symbolic Pauli exponent\n", + "\n", + "An arbitrary-angle Pauli exponent $e^{i\\alpha P} = \\cos(\\alpha)\\,I + i\\sin(\\alpha)\\,P$ is **not** a\n", + "stabilizer operation, so it cannot be applied as an ordinary Clifford gate. (Note that\n", + "`apply_pauli_exp` is only the *fixed* Clifford exponent $e^{i\\pi/4\\,P}$, not a free-angle one.)\n", + "\n", + "The simulator exposes it as a first-class operation parameterised by a **symbolic angle**: allocate the\n", + "angle with `allocate_symbolic_angle()`, then apply the exponent with `apply_symbolic_pauli_exp(P, angle)`.\n", + "The symbol $\\alpha$ stands for the angle of *every* $e^{i\\alpha P}$ at once. The simulator carries the\n", + "two halves of\n", + "$e^{i\\alpha P}|\\psi\\rangle = \\cos(\\alpha)\\,|\\psi\\rangle + i\\sin(\\alpha)\\,P|\\psi\\rangle$\n", + "together, with their **exact** relative phase, for all $\\alpha$ simultaneously. Two circuits agree for\n", + "every $\\alpha$ exactly when these halves match, phase included.\n", + "\n", + "Nothing here forms a $2^n$ state vector: phases are exact integer powers of $\\zeta_8 = e^{i\\pi/4}$ and\n", + "every step is polynomial in the number of qubits.\n", + "\n", + "> A symbolic angle is an **opaque handle**, not a number. Each one carries an `index` -- its\n", + "> subscript $k$ in $\\alpha_k$, fixed by allocation order -- and when two circuits are compared, angles\n", + "> sharing an `index` must correspond. Allocate a circuit's angles up front with\n", + "> `allocate_symbolic_angles(n)` and refer to them as `angles[k]`; the correspondence between two\n", + "> circuits is then explicit and independent of how either circuit is otherwise written." + ] + }, + { + "cell_type": "markdown", + "id": "221dcb60", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c2946815", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.715053Z", + "iopub.status.busy": "2026-07-05T18:33:59.714945Z", + "iopub.status.idle": "2026-07-05T18:33:59.718684Z", + "shell.execute_reply": "2026-07-05T18:33:59.717983Z" + } + }, + "outputs": [], + "source": [ + "from paulimer import (\n", + " PhasedOutcomeCompleteSimulation,\n", + " PhasedCircuitAction,\n", + " SparsePauli,\n", + " UnitaryOpcode,\n", + ")\n", + "\n", + "ZETA8_LABEL = {0: \"1\", 1: \"ζ₈\", 2: \"i\", 3: \"ζ₈³\", 4: \"-1\", 5: \"ζ₈⁵\", 6: \"-i\", 7: \"ζ₈⁷\"}" + ] + }, + { + "cell_type": "markdown", + "id": "06f95c97", + "metadata": {}, + "source": [ + "## A single symbolic Pauli exponent and its two halves\n", + "\n", + "We build $C_1\\, e^{i\\alpha Z_0}\\, C_2\\,|0\\rangle$ with $C_2 = H_0$ and $C_1 = H_0$. The simulator keeps\n", + "the identity half ($\\cos\\alpha$) and the $Z_0$ half ($i\\sin\\alpha$) of the exponent, each with its\n", + "exact $\\zeta_8$ phase." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "d6e65578", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.720189Z", + "iopub.status.busy": "2026-07-05T18:33:59.720112Z", + "iopub.status.idle": "2026-07-05T18:33:59.723242Z", + "shell.execute_reply": "2026-07-05T18:33:59.722332Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " identity half (weight cos a): zeta8^0 = 1\n", + " Z half (weight i*sin a): zeta8^0 = 1\n" + ] + } + ], + "source": [ + "sim = PhasedOutcomeCompleteSimulation(1)\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0]) # C_2\n", + "alpha = sim.allocate_symbolic_angle() # the exponent's symbolic angle alpha\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0\"), alpha) # e^{i alpha Z_0}\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0]) # C_1\n", + "\n", + "for fired in (0, 1):\n", + " exponent = sim.output_phase_exponent([bool(fired)])\n", + " half = \"identity half (weight cos a)\" if fired == 0 else \"Z half (weight i*sin a)\"\n", + " print(f\" {half}: zeta8^{exponent} = {ZETA8_LABEL[exponent]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "2e1ae37c", + "metadata": {}, + "source": [ + "## Verifying equivalence over all inputs with Choi states\n", + "\n", + "To check that two circuits implement the **same operator** on an *unknown* input, it is not enough to\n", + "run them on one fixed state: we must compare them on every input at once. Channel-state duality lets us\n", + "do this with a single stabilizer state. Entangle each of the $n$ system qubits with a fresh *reference*\n", + "qubit via a Bell pair $|\\Phi\\rangle$ (`prepare_bell_pairs` below), then apply the circuit to the system\n", + "qubits only. The resulting **Choi state** $(U \\otimes I)\\,|\\Phi\\rangle^{\\otimes n}$ determines $U$\n", + "completely.\n", + "\n", + "This is the same Bell-pair recipe used to verify Clifford circuits with `OutcomeCompleteSimulation`\n", + "(see *ZZ-measurement verification*); the only new ingredient is `apply_symbolic_pauli_exp`. Once the\n", + "circuit is applied, `sim.phased_action(input_qubits, output_qubits)` returns a `PhasedCircuitAction`,\n", + "and two actions are compared directly:\n", + "\n", + "- `a.is_equivalent(b)`: equal as operators on every input, **including** the exact relative phases\n", + " (up to one overall global phase);\n", + "- `a.is_equivalent_up_to_signs(b)`: equal only in their stabilizer (symplectic) action, ignoring all\n", + " phases, i.e. precisely what an ordinary, phase-blind stabilizer simulation sees.\n", + "\n", + "The comparison matches the symbolic angles by their **index**: $\\alpha_k$ of one circuit must\n", + "correspond to $\\alpha_k$ of the other." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a113dcd1", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.724997Z", + "iopub.status.busy": "2026-07-05T18:33:59.724828Z", + "iopub.status.idle": "2026-07-05T18:33:59.730661Z", + "shell.execute_reply": "2026-07-05T18:33:59.729163Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified: H . e^{i alpha Z} . H == e^{i alpha X} (as operators over all inputs)\n" + ] + } + ], + "source": [ + "def prepare_bell_pairs(sim, systems, references):\n", + " \"\"\"Entangle each system qubit with its reference qubit (the channel-state-duality setup).\"\"\"\n", + " for system, reference in zip(systems, references):\n", + " sim.apply_unitary(UnitaryOpcode.PrepareBell, [system, reference])\n", + "\n", + "\n", + "# H . e^{i alpha Z0} . H, as an operator on every input (system qubit 0, reference qubit 1).\n", + "sim = PhasedOutcomeCompleteSimulation(2)\n", + "prepare_bell_pairs(sim, [0], [1])\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0\"), alpha)\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0])\n", + "hzh = sim.phased_action([0], [0])\n", + "\n", + "# e^{i alpha X0}, on the same Choi layout.\n", + "sim = PhasedOutcomeCompleteSimulation(2)\n", + "prepare_bell_pairs(sim, [0], [1])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"X_0\"), alpha)\n", + "x_exp = sim.phased_action([0], [0])\n", + "\n", + "assert hzh.is_equivalent(x_exp)\n", + "print(\"verified: H . e^{i alpha Z} . H == e^{i alpha X} (as operators over all inputs)\")" + ] + }, + { + "cell_type": "markdown", + "id": "ef41d610", + "metadata": {}, + "source": [ + "## A multi-qubit, entangling equivalence\n", + "\n", + "The idiom extends verbatim to **entangling** Pauli exponents of arbitrary weight -- no new machinery.\n", + "A clean example: conjugating a single-qubit exponent by a `CNOT` turns it into a two-qubit\n", + "$ZZ$ exponent,\n", + "$$\\mathrm{CNOT}_{01}\\; e^{i\\alpha Z_1}\\; \\mathrm{CNOT}_{01} \\;=\\; e^{i\\alpha Z_0 Z_1},$$\n", + "because $\\mathrm{CNOT}_{01}$ conjugates $Z_1 \\mapsto Z_0 Z_1$. We verify it as operators (Choi state,\n", + "$n=2$), and confirm that dropping the conjugation -- a bare $e^{i\\alpha Z_1}$ -- is genuinely different." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "3cbf62d9", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.733195Z", + "iopub.status.busy": "2026-07-05T18:33:59.733005Z", + "iopub.status.idle": "2026-07-05T18:33:59.738664Z", + "shell.execute_reply": "2026-07-05T18:33:59.738202Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified: e^{i alpha Z0Z1} == CNOT01 . e^{i alpha Z1} . CNOT01 (entangling)\n", + "verified: e^{i alpha Z0Z1} != e^{i alpha Z1} (the CNOT conjugation genuinely matters)\n" + ] + } + ], + "source": [ + "# e^{i alpha Z0 Z1} (system qubits 0,1; reference qubits 2,3)\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0 Z_1\"), alpha)\n", + "zz_direct = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "# CNOT01 . e^{i alpha Z1} . CNOT01\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_1\"), alpha)\n", + "sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1])\n", + "zz_via_cnot = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "# e^{i alpha Z1} alone (conjugation dropped)\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_1\"), alpha)\n", + "z1_only = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "assert zz_direct.is_equivalent(zz_via_cnot)\n", + "print(\"verified: e^{i alpha Z0Z1} == CNOT01 . e^{i alpha Z1} . CNOT01 (entangling)\")\n", + "\n", + "assert not zz_direct.is_equivalent(z1_only)\n", + "print(\"verified: e^{i alpha Z0Z1} != e^{i alpha Z1} (the CNOT conjugation genuinely matters)\")" + ] + }, + { + "cell_type": "markdown", + "id": "e11fb15c", + "metadata": {}, + "source": [ + "## A mixed-basis Pauli exponent\n", + "\n", + "Symbolic exponents are not restricted to the $Z$ basis: the observable $P$ can be **any** Pauli.\n", + "Conjugating the $ZZ$ exponent above by a Hadamard on qubit $0$ rotates $Z_0 \\mapsto X_0$, giving a\n", + "mixed-basis, weight-two exponent\n", + "$$H_0\\; e^{i\\alpha Z_0 Z_1}\\; H_0 \\;=\\; e^{i\\alpha X_0 Z_1}.$$\n", + "We verify this equivalence as operators, and confirm that the un-conjugated $e^{i\\alpha Z_0 Z_1}$ is a\n", + "genuinely different channel -- the choice of basis in the exponent matters." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "c0063c8e", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.740797Z", + "iopub.status.busy": "2026-07-05T18:33:59.740605Z", + "iopub.status.idle": "2026-07-05T18:33:59.746312Z", + "shell.execute_reply": "2026-07-05T18:33:59.745161Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified: e^{i alpha X0 Z1} == H0 . e^{i alpha Z0 Z1} . H0 (mixed-basis, higher weight)\n", + "verified: e^{i alpha X0 Z1} != e^{i alpha Z0 Z1} (the exponent's basis matters)\n" + ] + } + ], + "source": [ + "# e^{i alpha X0 Z1} -- a mixed-basis (non-Z), higher-weight Pauli exponent\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"X_0 Z_1\"), alpha)\n", + "xz_direct = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "# H0 . e^{i alpha Z0 Z1} . H0\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0 Z_1\"), alpha)\n", + "sim.apply_unitary(UnitaryOpcode.Hadamard, [0])\n", + "xz_via_h = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "# the un-conjugated ZZ exponent is a different channel\n", + "sim = PhasedOutcomeCompleteSimulation(4)\n", + "prepare_bell_pairs(sim, [0, 1], [2, 3])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0 Z_1\"), alpha)\n", + "zz_only = sim.phased_action([0, 1], [0, 1])\n", + "\n", + "assert xz_direct.is_equivalent(xz_via_h)\n", + "print(\"verified: e^{i alpha X0 Z1} == H0 . e^{i alpha Z0 Z1} . H0 (mixed-basis, higher weight)\")\n", + "\n", + "assert not xz_direct.is_equivalent(zz_only)\n", + "print(\"verified: e^{i alpha X0 Z1} != e^{i alpha Z0 Z1} (the exponent's basis matters)\")" + ] + }, + { + "cell_type": "markdown", + "id": "aa7233d2", + "metadata": {}, + "source": [ + "## A phase difference that ordinary simulation misses\n", + "\n", + "Finally, a difference that is *purely* a phase. The exponents $e^{+i\\alpha Z}$ and $e^{-i\\alpha Z}$\n", + "differ only by the sign of the Pauli, $+Z$ versus $-Z$. Since $+Z$ and $-Z$ have the **same symplectic\n", + "action**, an ordinary stabilizer simulation cannot tell them apart -- *even* under the Choi comparison\n", + "above. Yet they are physically different, differing by a relative $-1$ on the $Z$ half of the exponent.\n", + "The two comparisons make this explicit: `is_equivalent_up_to_signs` (phase-blind) reports them equal,\n", + "while `is_equivalent` (phase-aware) separates them." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "61093f71", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.748049Z", + "iopub.status.busy": "2026-07-05T18:33:59.747876Z", + "iopub.status.idle": "2026-07-05T18:33:59.752269Z", + "shell.execute_reply": "2026-07-05T18:33:59.751210Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ignoring phase (is_equivalent_up_to_signs): e^{+i alpha Z} and e^{-i alpha Z} INDISTINGUISHABLE\n", + "Tracking phase (is_equivalent): e^{+i alpha Z} != e^{-i alpha Z}\n" + ] + } + ], + "source": [ + "# e^{+i alpha Z0}\n", + "sim = PhasedOutcomeCompleteSimulation(2)\n", + "prepare_bell_pairs(sim, [0], [1])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"Z_0\"), alpha)\n", + "positive = sim.phased_action([0], [0])\n", + "\n", + "# e^{-i alpha Z0}\n", + "sim = PhasedOutcomeCompleteSimulation(2)\n", + "prepare_bell_pairs(sim, [0], [1])\n", + "alpha = sim.allocate_symbolic_angle()\n", + "sim.apply_symbolic_pauli_exp(SparsePauli(\"-Z_0\"), alpha)\n", + "negative = sim.phased_action([0], [0])\n", + "\n", + "assert positive.is_equivalent_up_to_signs(negative)\n", + "print(\"Ignoring phase (is_equivalent_up_to_signs): e^{+i alpha Z} and e^{-i alpha Z} INDISTINGUISHABLE\")\n", + "\n", + "assert not positive.is_equivalent(negative)\n", + "print(\"Tracking phase (is_equivalent): e^{+i alpha Z} != e^{-i alpha Z}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "36d37f2f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.754274Z", + "iopub.status.busy": "2026-07-05T18:33:59.754168Z", + "iopub.status.idle": "2026-07-05T18:33:59.756894Z", + "shell.execute_reply": "2026-07-05T18:33:59.756383Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "e^{+i alpha Z} half phases: exponents (0, 0) = ('1', '1')\n", + "e^{-i alpha Z} half phases: exponents (0, 4) = ('1', '-1')\n", + "\n", + "The Z half differs by zeta8^4 = -1, exactly the relative phase ordinary simulation discards.\n" + ] + } + ], + "source": [ + "def exponent_half_phases(observable):\n", + " \"\"\"The exact zeta8 exponent on each half (identity, P) of e^{i alpha observable}.\"\"\"\n", + " sim = PhasedOutcomeCompleteSimulation(2)\n", + " prepare_bell_pairs(sim, [0], [1])\n", + " alpha = sim.allocate_symbolic_angle()\n", + " sim.apply_symbolic_pauli_exp(observable, alpha)\n", + " return tuple(sim.output_phase_exponent([bool(fired)]) for fired in (0, 1))\n", + "\n", + "\n", + "for name, observable in ((\"e^{+i alpha Z}\", SparsePauli(\"Z_0\")), (\"e^{-i alpha Z}\", SparsePauli(\"-Z_0\"))):\n", + " phases = exponent_half_phases(observable)\n", + " labels = tuple(ZETA8_LABEL[power] for power in phases)\n", + " print(f\"{name} half phases: exponents {phases} = {labels}\")\n", + "print(\"\\nThe Z half differs by zeta8^4 = -1, exactly the relative phase ordinary simulation discards.\")" + ] + }, + { + "cell_type": "markdown", + "id": "663bbabe", + "metadata": {}, + "source": [ + "## Executing a diagonal channel remotely: \"ejection\"\n", + "\n", + "A channel that is **diagonal in the $Z$ basis** -- any number of $Z$-Pauli exponents together with\n", + "non-destructive $Z$-parity measurements -- can be executed **remotely**, the way a measurement-based or\n", + "lattice-surgery gadget would. Copy each system qubit onto a fresh ancilla with a `CNOT`, run the\n", + "diagonal channel on the **ancillas**, then **measure** each ancilla in the $X$ basis; a $-$ outcome\n", + "triggers a conditional $Z$ correction on the corresponding system qubit. Averaging over the measurement\n", + "outcomes, the channel the gadget implements on the inputs is exactly the direct channel.\n", + "\n", + "The example below ejects a genuinely non-trivial three-qubit channel: three overlapping $Z$-Pauli\n", + "exponents,\n", + "$$e^{i\\alpha Z_0}, \\qquad e^{i\\beta Z_1 Z_2}, \\qquad e^{i\\gamma Z_0 Z_1},$$\n", + "followed by a non-destructive measurement of the three-qubit parity $Z_0 Z_1 Z_2$. This mixes three\n", + "kinds of bit at once: the symbolic angles of the exponents (matched one-to-one against the direct\n", + "circuit), the non-destructive parity outcome (a measurement outcome, observed identically in both circuits), and\n", + "the destructive $X$ read-outs that drive the corrections (measurement outcomes, marginalized away). The phased\n", + "comparison certifies the ejected gadget equals the direct channel for all $\\alpha, \\beta, \\gamma$." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "0ed537fd", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T18:33:59.758657Z", + "iopub.status.busy": "2026-07-05T18:33:59.758563Z", + "iopub.status.idle": "2026-07-05T18:33:59.763579Z", + "shell.execute_reply": "2026-07-05T18:33:59.762789Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "verified: three-qubit Z-diagonal channel (3 Pauli exponents + a 3-qubit parity measurement)\n", + " ejected through ancillas == the channel applied directly\n" + ] + } + ], + "source": [ + "def apply_z_diagonal_channel(sim, support, angles):\n", + " \"\"\"A Z-diagonal channel on three qubits: three overlapping Pauli exponents e^{i alpha_k P_k}\n", + " parameterised by the given symbolic angles, then a non-destructive measurement of the three-qubit\n", + " Z parity. The k-th exponent uses angles[k], so two circuits given angles with matching index apply\n", + " the same channel.\"\"\"\n", + " q0, q1, q2 = support\n", + " observables = (\n", + " SparsePauli(f\"Z_{q0}\"), # angles[0] . Z_q0\n", + " SparsePauli(f\"Z_{q1} Z_{q2}\"), # angles[1] . Z_q1 Z_q2\n", + " SparsePauli(f\"Z_{q0} Z_{q1}\"), # angles[2] . Z_q0 Z_q1\n", + " )\n", + " for angle, observable in zip(angles, observables):\n", + " sim.apply_symbolic_pauli_exp(observable, angle)\n", + " sim.measure(SparsePauli(f\"Z_{q0} Z_{q1} Z_{q2}\")) # non-destructive 3-qubit parity measurement\n", + "\n", + "\n", + "system = [0, 1, 2] # system qubits (their references are 3, 4, 5)\n", + "ancilla = [6, 7, 8]\n", + "\n", + "# Direct: run the diagonal channel straight on the system qubits.\n", + "direct = PhasedOutcomeCompleteSimulation(6)\n", + "prepare_bell_pairs(direct, system, [3, 4, 5])\n", + "angles = direct.allocate_symbolic_angles(3) # alpha_0, alpha_1, alpha_2\n", + "apply_z_diagonal_channel(direct, system, angles)\n", + "direct_action = direct.phased_action(system, system)\n", + "\n", + "# Ejected: copy each system qubit onto an ancilla, run the channel on the ancillas, then read the\n", + "# ancillas out in X and correct. Allocating the angles the same way makes alpha_k here the same alpha_k\n", + "# as in the direct circuit -- the comparison pairs them by index.\n", + "ejected = PhasedOutcomeCompleteSimulation(9)\n", + "prepare_bell_pairs(ejected, system, [3, 4, 5])\n", + "angles = ejected.allocate_symbolic_angles(3) # the same alpha_0, alpha_1, alpha_2\n", + "for system_qubit, ancilla_qubit in zip(system, ancilla):\n", + " ejected.apply_unitary(UnitaryOpcode.ControlledX, [system_qubit, ancilla_qubit])\n", + "apply_z_diagonal_channel(ejected, ancilla, angles)\n", + "for system_qubit, ancilla_qubit in zip(system, ancilla):\n", + " outcome = ejected.measure(SparsePauli(f\"X_{ancilla_qubit}\"))\n", + " ejected.apply_conditional_pauli(SparsePauli(f\"Z_{system_qubit}\"), [outcome])\n", + "ejected_action = ejected.phased_action(system, system)\n", + "\n", + "assert direct_action.is_equivalent(ejected_action)\n", + "print(\"verified: three-qubit Z-diagonal channel (3 Pauli exponents + a 3-qubit parity measurement)\")\n", + "print(\" ejected through ancillas == the channel applied directly\")" + ] + }, + { + "cell_type": "markdown", + "id": "c4b50018", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "- A symbolic Pauli exponent $e^{i\\alpha P}$ is added by allocating an opaque angle handle with\n", + " `allocate_symbolic_angle()` (or a batch with `allocate_symbolic_angles(n)`) and applying\n", + " `apply_symbolic_pauli_exp(P, angle)`. This works for an **arbitrary** Pauli $P$ of any weight:\n", + " multi-qubit and entangling exponents need nothing new, and a single angle may parameterise several\n", + " exponents to model a shared $\\alpha$. Two circuits' angles correspond when they share an `index`\n", + " ($\\alpha_k$), so allocating each circuit's angles up front and indexing them keeps the comparison\n", + " explicit.\n", + "- To compare two circuits as **operators** on an unknown input, build their **Choi states**\n", + " (`prepare_bell_pairs` to entangle every system qubit with a reference, then apply the circuit to the\n", + " system qubits) and call `sim.phased_action(...)`. This is the same Bell-pair recipe used to verify\n", + " Clifford circuits with `OutcomeCompleteSimulation`. The resulting `PhasedCircuitAction` objects\n", + " compare with `is_equivalent` (phase-aware) and `is_equivalent_up_to_signs` (phase-blind).\n", + "- The phased simulator tracks both halves of every exponent with their **exact** $\\zeta_8$ phase,\n", + " precisely the information ordinary stabilizer simulation throws away. Here it is the only thing\n", + " distinguishing $e^{+i\\alpha Z}$ from $e^{-i\\alpha Z}$, whose Paulis $+Z$ and $-Z$ share a symplectic\n", + " action.\n", + "- A whole $Z$-diagonal channel -- several Pauli exponents plus non-destructive parity measurements --\n", + " can be **ejected** onto ancillas and realised through $X$ read-outs and Pauli corrections; the phased\n", + " comparison certifies the resulting channel equals the direct one for every choice of angles." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/paulimer/bindings/python/paulimer.pyi b/paulimer/bindings/python/paulimer.pyi index ee79347d..bc4d4c85 100644 --- a/paulimer/bindings/python/paulimer.pyi +++ b/paulimer/bindings/python/paulimer.pyi @@ -27,7 +27,10 @@ __all__ = [ "PauliDistribution", "PauliFault", "PauliGroup", + "PhasedCircuitAction", + "PhasedOutcomeCompleteSimulation", "SparsePauli", + "SymbolicAngle", "UnitaryOpcode", "centralizer_of", "encoding_clifford_of", @@ -644,6 +647,16 @@ class CliffordUnitary: """Split into phased CSS components.""" ... + def to_pauli_exponents(self) -> list[SparsePauli]: + """Decompose into an ordered product of pi/4 Pauli exponents. + + Returns Paulis ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, ..., + ``exp(i pi/4 P_k)`` to the identity via ``left_mul_pauli_exp`` reproduces this Clifford + exactly, including the Pauli-image signs. The sign of each returned Pauli selects + ``exp(+i pi/4 P)`` or ``exp(-i pi/4 P)``. + """ + ... + def __mul__(self, other: "CliffordUnitary", /) -> "CliffordUnitary": ... def left_mul(self, unitary_op: UnitaryOpcode, support: Sequence[int]) -> None: ... def left_mul_clifford( @@ -1030,6 +1043,306 @@ class OutcomeCompleteSimulation: """ ... +@final +class SymbolicAngle: + """An opaque handle to a symbolic angle ``alpha`` of a parameterised circuit. + + A symbolic angle is the free parameter of a Pauli exponent ``e^{i alpha P}``. Obtain one from + :meth:`PhasedOutcomeCompleteSimulation.allocate_symbolic_angle` (or a batch from + :meth:`PhasedOutcomeCompleteSimulation.allocate_symbolic_angles`) and pass it to + :meth:`PhasedOutcomeCompleteSimulation.apply_symbolic_pauli_exp`. The handle is opaque: its only + observable feature is its :attr:`index`, the angle's subscript ``k`` in ``alpha_k``, fixed by the + order in which angles are allocated. When two circuits are compared with + :meth:`PhasedOutcomeCompleteSimulation.phased_action`, angles with the same index are required to + correspond, so describing both circuits in terms of the ``k``-th angle is what makes the + comparison meaningful -- regardless of how the rest of each circuit is written. + """ + + @property + def index(self) -> int: + """The subscript ``k`` identifying this angle as ``alpha_k``, set by allocation order.""" + ... + + def __eq__(self, other: object, /) -> bool: ... + def __hash__(self) -> int: ... + def __repr__(self) -> str: ... + +@final +class PhasedOutcomeCompleteSimulation: + """Outcome-complete stabilizer simulation that also tracks the exact global phase. + + This is the global-phase-resolving generalization of + :class:`OutcomeCompleteSimulation`, implementing Algorithm 4.2 of + arXiv:2603.24717 ("phased outcome-complete simulation"). Like its phaseless + counterpart it tracks all ``2^n_random`` measurement branches simultaneously, + but the encoded state is maintained with its *exact* global phase rather than + only up to a global phase. This enables exact equality checking of non-stabilizer + circuits (e.g. circuits with symbolic single-qubit Pauli exponents). + + For a random-bit assignment ``r`` the encoded state is + + ``i^ (-1)^ R|A r>`` + + where ``R`` is the phased state encoder, ``A`` the sign matrix, ``B`` the quadratic + phase matrix, ``p`` the linear ``i``-phase vector, and ``s`` the linear ``-1``-phase + vector. The scalar prefactor is exposed via :meth:`output_phase_exponent`. + + Examples: + >>> sim = PhasedOutcomeCompleteSimulation(2) + >>> sim.apply_unitary(UnitaryOpcode.Hadamard, [0]) + >>> sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + >>> sim.measure(SparsePauli("X_0")) + >>> exponent = sim.output_phase_exponent([True]) # zeta_8 exponent for r = (1,) + """ + + def __new__(cls, qubit_count: int = 0) -> "PhasedOutcomeCompleteSimulation": + """Create a simulation with the specified number of qubits.""" + ... + + @property + def qubit_count(self) -> int: ... + @property + def qubit_capacity(self) -> int: ... + @property + def outcome_count(self) -> int: ... + @property + def outcome_capacity(self) -> int: ... + @property + def random_outcome_count(self) -> int: ... + @property + def random_outcome_capacity(self) -> int: ... + @property + def random_bit_count(self) -> int: ... + def apply_unitary( + self, unitary_op: UnitaryOpcode, support: Sequence[int] + ) -> None: ... + def apply_pauli_exp(self, observable: SparsePauli) -> None: ... + def apply_pauli( + self, observable: SparsePauli, controlled_by: SparsePauli | None = None + ) -> None: ... + def apply_conditional_pauli( + self, + observable: SparsePauli, + outcomes: Sequence[int], + parity: bool = True, + ) -> None: ... + def apply_permutation( + self, permutation: Sequence[int], supported_by: Sequence[int] | None = None + ) -> None: ... + def apply_clifford( + self, clifford: CliffordUnitary, supported_by: Sequence[int] | None = None + ) -> None: + """Unsupported on the phased simulator. + + A phaseless :class:`CliffordUnitary` does not determine the exact global phase that this + simulator tracks. Apply Cliffords through :meth:`apply_unitary`, :meth:`apply_pauli`, or + :meth:`apply_pauli_exp` instead. + + Raises: + NotImplementedError: Always. + """ + ... + def measure( + self, observable: SparsePauli, hint: SparsePauli | None = None + ) -> int: ... + def allocate_random_bit(self) -> int: ... + def reserve_qubits(self, new_qubit_capacity: int) -> None: ... + def reserve_outcomes( + self, new_outcome_capacity: int, new_random_outcome_capacity: int + ) -> None: ... + def is_stabilizer( + self, + observable: SparsePauli, + ignore_sign: bool = False, + sign_parity: Sequence[int] = ..., # type: ignore[assignment] + ) -> bool: + """Check if an observable is a stabilizer of the current state.""" + ... + + @staticmethod + def with_capacity( + num_qubits: int, num_outcomes: int, num_random_outcomes: int + ) -> "PhasedOutcomeCompleteSimulation": + """Create simulation with pre-allocated capacity.""" + ... + + @property + def random_outcome_indicator(self) -> BitVector: + """Indicator of which outcomes are random (vs deterministic).""" + ... + + @property + def clifford(self) -> CliffordUnitary: + """Clifford unitary encoding the current stabilizer state (global phase discarded).""" + ... + + @property + def sign_matrix(self) -> BitMatrix: + """Sign matrix A mapping random outcomes to the computational-basis register. + + Shape: (qubit_count, random_outcome_count) + """ + ... + + @property + def quadratic_phase_matrix(self) -> BitMatrix: + """Quadratic phase matrix B contributing the (-1)^ factor. + + Shape: (random_outcome_count, random_outcome_count) + """ + ... + + @property + def outcome_matrix(self) -> BitMatrix: + """Outcome matrix M encoding all 2^k measurement branches. + + Shape: (outcome_count, random_outcome_count) + """ + ... + + @property + def outcome_shift(self) -> BitVector: + """Outcome shift vector v_0 representing deterministic outcome contributions. + + Length: outcome_count + """ + ... + + @property + def linear_i_phase(self) -> BitVector: + """Linear i-phase vector p contributing the i^ factor. + + Length: random_outcome_count + """ + ... + + @property + def linear_sign_phase(self) -> BitVector: + """Linear sign-phase vector s contributing the (-1)^ factor. + + Length: random_outcome_count + """ + ... + + def output_phase_exponent(self, random_bits: Sequence[bool]) -> int: + """Return the zeta_8 = e^{i pi/4} exponent of the scalar prefactor. + + For the given random-bit assignment ``r`` this is the exponent (modulo 8) of + the scalar ``i^ (-1)^`` multiplying ``R|A r>`` in the output + state. The phase of ``R|A r>`` itself is carried by the phased encoder. + + Args: + random_bits: Boolean assignment for each random outcome (length at least + ``random_outcome_count``). + """ + ... + + def allocate_symbolic_angle(self) -> SymbolicAngle: + """Allocate a fresh symbolic angle ``alpha``. + + Returns an opaque :class:`SymbolicAngle` handle; pass it to + :meth:`apply_symbolic_pauli_exp` to apply ``e^{i alpha P}``. Angles are numbered by + allocation order (the handle's :attr:`SymbolicAngle.index`), and when two circuits are + compared with :meth:`phased_action` the angle with a given index in one must correspond to + the same index in the other. To allocate several at once, use + :meth:`allocate_symbolic_angles`. + """ + ... + + def allocate_symbolic_angles(self, count: int) -> list[SymbolicAngle]: + """Allocate ``count`` fresh symbolic angles ``alpha_0, ..., alpha_{count-1}`` at once. + + Returns the :class:`SymbolicAngle` handles in order, so ``angles[k]`` is ``alpha_k``. + Allocating all of a circuit's angles up front and referring to them by index keeps the + correspondence between two circuits explicit and independent of how either is otherwise + written. + """ + ... + + @property + def symbolic_angles(self) -> list[SymbolicAngle]: + """All symbolic angles allocated so far, in order (``angles[k]`` is ``alpha_k``).""" + ... + + def apply_symbolic_pauli_exp(self, observable: SparsePauli, angle: SymbolicAngle) -> None: + """Apply a symbolic Pauli exponent ``e^{i alpha P}`` parameterised by ``angle``. + + ``angle`` must be a :class:`SymbolicAngle` obtained from :meth:`allocate_symbolic_angle` + or :meth:`allocate_symbolic_angles`. This is the high-level way to add a free-angle + exponent ``e^{i alpha P}`` for an arbitrary Pauli ``P``. The same ``angle`` may parameterise + several exponents (a shared ``alpha``), and angles with matching index in two circuits are + what make those circuits' exponents correspond when their phased actions are compared. + """ + ... + + def phased_action( + self, input_qubits: Sequence[int], output_qubits: Sequence[int] + ) -> PhasedCircuitAction: + """Compute the phased Choi action of the circuit recorded in this simulation. + + Returns a :class:`PhasedCircuitAction` capturing how the circuit acts on every + input at once, including the exact relative phases between measurement branches. + This is the global-phase-resolving counterpart of the (phaseless) circuit action + used to compare stabilizer circuits. + + Before calling this, the Choi state must already be prepared: entangle each + ``input_qubits[k]`` with a fresh reference qubit via + ``UnitaryOpcode.PrepareBell`` and then apply the circuit to the system qubits + only. The reference qubit for ``input_qubits[k]`` is ``system_qubit_count + k``, + where ``system_qubit_count`` is one past the largest index in ``input_qubits`` or + ``output_qubits`` (so for ``n`` system qubits ``0..n`` the references are + ``n..2n``). + + Symbolic angles (allocated with :meth:`allocate_symbolic_angle`) are matched + one-to-one by index between the two compared actions, while genuine measurement + randomness is marginalized over (see :meth:`PhasedCircuitAction.is_equivalent`). + + Args: + input_qubits: System qubits entangled with reference qubits. + output_qubits: System qubits carrying the circuit's output. + + Raises: + ValueError: If the non-output system qubits remain entangled with the rest of + the state. + """ + ... + +@final +class PhasedCircuitAction: + """The action of a circuit on every input, with exact relative branch phases. + + Produced by :meth:`PhasedOutcomeCompleteSimulation.phased_action`. Two actions are + compared up to a single overall global phase; the *relative* phases between branches + are retained, so circuits that act identically on the Pauli group but differ by a + branch-dependent phase (for example ``e^{i a Z}`` versus ``e^{-i a Z}``, whose + conditioned Paulis ``+Z`` and ``-Z`` share a symplectic action) are distinguished. + + Symbolic angles are matched one-to-one by index between the two compared actions, while + genuine measurement random bits are marginalized over (as in the phaseless comparison). + """ + + @property + def choi_state_stabilizers(self) -> list[SparsePauli]: + """Canonical stabilizers of the circuit's Choi state.""" + ... + + def is_equivalent(self, other: PhasedCircuitAction) -> bool: + """Whether two circuits implement the same operator on every input. + + Compares both the stabilizer (symplectic) action and the exact relative branch + phases, up to a single global phase. Symbolic angles are matched one-to-one by + index with those of ``other``, while genuine measurement randomness is marginalized. + """ + ... + + def is_equivalent_up_to_signs(self, other: PhasedCircuitAction) -> bool: + """Whether two circuits agree on their stabilizer action, ignoring all phases. + + This is the phaseless comparison; use :meth:`is_equivalent` to additionally + require the relative branch phases to match. + """ + ... + @final class OutcomeFreeSimulation: """Stabilizer simulation without tracking specific measurement outcomes. diff --git a/paulimer/bindings/python/src/lib.rs b/paulimer/bindings/python/src/lib.rs index 394a0723..2053ad9f 100644 --- a/paulimer/bindings/python/src/lib.rs +++ b/paulimer/bindings/python/src/lib.rs @@ -22,7 +22,10 @@ pub use py_frame_propagator::PyFramePropagator; pub use py_noise::{PyFault, PyOutcomeCondition, PyPauliDistribution}; pub use py_pauli_group::{py_centralizer_of, py_symplectic_form_of, PyPauliGroup}; pub use py_sparse_pauli::PySparsePauli; -pub use simulation::{PyOutcomeCompleteSimulation, PyOutcomeFreeSimulation, PyOutcomeSpecificSimulation}; +pub use simulation::{ + PyOutcomeCompleteSimulation, PyOutcomeFreeSimulation, PyOutcomeSpecificSimulation, PyPhasedCircuitAction, + PyPhasedOutcomeCompleteSimulation, PySymbolicAngle, +}; /// # Errors /// @@ -35,6 +38,9 @@ pub fn paulimer(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; + m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/paulimer/bindings/python/src/py_clifford.rs b/paulimer/bindings/python/src/py_clifford.rs index bdf51d2a..01a026f3 100644 --- a/paulimer/bindings/python/src/py_clifford.rs +++ b/paulimer/bindings/python/src/py_clifford.rs @@ -1,7 +1,7 @@ use derive_more::{Deref, DerefMut, From, Into}; use paulimer::clifford::{ - group_encoding_clifford_of, split_phased_css, split_qubit_cliffords_and_css, Clifford, CliffordMutable, - CliffordUnitary, XOrZ, + clifford_to_pauli_exponents, group_encoding_clifford_of, split_phased_css, split_qubit_cliffords_and_css, Clifford, + CliffordMutable, CliffordUnitary, XOrZ, }; use paulimer::pauli::{as_sparse, DensePauli, SparsePauli}; use pyo3::exceptions::PyValueError; @@ -162,6 +162,13 @@ impl PyCliffordUnitary { Some((left.into(), right.into())) } + fn to_pauli_exponents(&self) -> Vec { + clifford_to_pauli_exponents(&self.inner) + .into_iter() + .map(PySparsePauli::from) + .collect() + } + fn preimage_x(&self, qubit_index: usize) -> PyDensePauli { PyDensePauli { inner: self.inner.preimage_x(qubit_index), diff --git a/paulimer/bindings/python/src/simulation.rs b/paulimer/bindings/python/src/simulation.rs index 8e00d408..69a65b10 100644 --- a/paulimer/bindings/python/src/simulation.rs +++ b/paulimer/bindings/python/src/simulation.rs @@ -2,10 +2,13 @@ use std::ops::{Deref, DerefMut}; use binar::{BitMatrix, BitVec}; use paulimer::clifford::CliffordUnitary; +use pauliverse::action::{phased_action_from_simulation, PhasedCircuitAction}; use pauliverse::outcome_complete_simulation::OutcomeCompleteSimulation; use pauliverse::outcome_free_simulation::OutcomeFreeSimulation; use pauliverse::outcome_specific_simulation::OutcomeSpecificSimulation; +use pauliverse::phased_outcome_complete_simulation::PhasedOutcomeCompleteSimulation; use pauliverse::Simulation; +use pyo3::exceptions::{PyNotImplementedError, PyValueError}; use pyo3::prelude::*; use crate::enums::PyUnitaryOp; @@ -33,8 +36,57 @@ pub struct PyOutcomeFreeSimulation { inner: OutcomeFreeSimulation, } +#[derive(derive_more::Deref, derive_more::DerefMut, derive_more::From)] +#[must_use] +#[pyclass(name = "PhasedOutcomeCompleteSimulation", module = "paulimer")] +pub struct PyPhasedOutcomeCompleteSimulation { + inner: PhasedOutcomeCompleteSimulation, +} + +/// An opaque handle to a symbolic angle `α` of a parameterised circuit. +/// +/// A symbolic angle is the free parameter of a Pauli exponent `e^{iα P}`. Obtain one from +/// [`PhasedOutcomeCompleteSimulation.allocate_symbolic_angle`] (or a batch from +/// `allocate_symbolic_angles`) and pass it to `apply_symbolic_pauli_exp`. The handle is opaque: +/// its only observable feature is its [`index`], the angle's subscript `k` in `α_k`, fixed by the +/// order in which angles are allocated. When two circuits are compared with `phased_action`, angles +/// with the same `index` are required to correspond, so describing both circuits in terms of the +/// `k`-th angle is what makes the comparison meaningful -- regardless of how the rest of each +/// circuit is written. +#[derive(Clone)] +#[pyclass(name = "SymbolicAngle", module = "paulimer", frozen)] +pub struct PySymbolicAngle { + outcome: usize, + index: usize, +} + +#[pymethods] +impl PySymbolicAngle { + /// The subscript `k` identifying this angle as `α_k`, set by allocation order. + #[getter] + #[must_use] + pub fn index(&self) -> usize { + self.index + } + + #[must_use] + pub fn __repr__(&self) -> String { + format!("SymbolicAngle(index={})", self.index) + } + + #[must_use] + pub fn __eq__(&self, other: &Self) -> bool { + self.outcome == other.outcome + } + + #[must_use] + pub fn __hash__(&self) -> u64 { + self.outcome as u64 + } +} + macro_rules! impl_simulation { - ($struct_name:ty, $wrapper_struct:ty { $($inside:tt)* }) => { + ($struct_name:ty, $wrapper_struct:ty, clifford_supported = $clifford_supported:literal { $($inside:tt)* }) => { #[pymethods] impl $wrapper_struct { #[new] @@ -116,11 +168,19 @@ macro_rules! impl_simulation { Simulation::permute(self.deref_mut(), &permutation, &support); } - #[allow(clippy::needless_pass_by_value)] + #[allow(clippy::needless_pass_by_value, clippy::missing_errors_doc)] #[pyo3(signature=(clifford, supported_by=None))] - pub fn apply_clifford(&mut self, clifford: &PyCliffordUnitary, supported_by: Option>) { + pub fn apply_clifford(&mut self, clifford: &PyCliffordUnitary, supported_by: Option>) -> PyResult<()> { + if !$clifford_supported { + return Err(PyNotImplementedError::new_err( + "this simulator tracks the exact global phase, which a phaseless CliffordUnitary \ + does not determine; apply Cliffords through apply_unitary, apply_pauli, or \ + apply_pauli_exp instead", + )); + } let support = supported_by.unwrap_or_else(|| (0..self.deref().qubit_count()).collect()); Simulation::clifford(self.deref_mut(), &clifford.inner, &support); + Ok(()) } #[pyo3(signature=(observable, hint=None))] @@ -177,7 +237,8 @@ macro_rules! impl_simulation { impl_simulation!( OutcomeCompleteSimulation, - PyOutcomeCompleteSimulation { + PyOutcomeCompleteSimulation, + clifford_supported = true { #[getter] pub fn clifford(&self) -> PyCliffordUnitary { PyCliffordUnitary { @@ -203,7 +264,8 @@ impl_simulation!( impl_simulation!( OutcomeFreeSimulation, - PyOutcomeFreeSimulation { + PyOutcomeFreeSimulation, + clifford_supported = true { #[getter] pub fn clifford(&self) -> PyCliffordUnitary { let c: CliffordUnitary = self.deref().state_encoder().clone().into(); @@ -213,7 +275,8 @@ impl_simulation!( impl_simulation!( OutcomeSpecificSimulation, - PyOutcomeSpecificSimulation { + PyOutcomeSpecificSimulation, + clifford_supported = true { #[getter] pub fn clifford(&self) -> PyCliffordUnitary { PyCliffordUnitary { @@ -238,3 +301,143 @@ impl_simulation!( OutcomeSpecificSimulation::new_with_seeded_random_outcomes(num_qubits, seed).into() } }); + +impl_simulation!( + PhasedOutcomeCompleteSimulation, + PyPhasedOutcomeCompleteSimulation, + clifford_supported = false { + #[getter] + pub fn clifford(&self) -> PyCliffordUnitary { + PyCliffordUnitary { + inner: self.deref().state_encoder(), + } + } + + #[getter] + pub fn sign_matrix(&self) -> BitMatrix { + self.inner.sign_matrix() + } + + #[getter] + pub fn quadratic_phase_matrix(&self) -> BitMatrix { + self.inner.quadratic_phase_matrix() + } + + #[getter] + pub fn outcome_matrix(&self) -> BitMatrix { + self.inner.outcome_matrix() + } + + #[getter] + pub fn outcome_shift(&self) -> BitVec { + self.inner.outcome_shift() + } + + #[getter] + pub fn linear_i_phase(&self) -> BitVec { + self.inner.linear_i_phase() + } + + #[getter] + pub fn linear_sign_phase(&self) -> BitVec { + self.inner.linear_sign_phase() + } + + #[allow(clippy::needless_pass_by_value)] + #[must_use] + pub fn output_phase_exponent(&self, random_bits: Vec) -> u8 { + self.inner.output_phase_exponent(&random_bits) + } + + /// Allocate a fresh symbolic angle `α`. + /// + /// Returns an opaque [`SymbolicAngle`] handle; pass it to [`apply_symbolic_pauli_exp`] to + /// apply `e^{iα P}`. Angles are numbered by allocation order (the returned handle's + /// `index`), and when two circuits are compared with `phased_action` the angle with a given + /// `index` in one must correspond to the same `index` in the other. To allocate several at + /// once, use [`allocate_symbolic_angles`]. + pub fn allocate_symbolic_angle(&mut self) -> PySymbolicAngle { + let index = self.inner.symbolic_angle_indicator().iter().filter(|&&is_angle| is_angle).count(); + let outcome = self.inner.allocate_symbolic_angle(); + PySymbolicAngle { outcome, index } + } + + /// Allocate `count` fresh symbolic angles `α_0, ..., α_{count-1}` at once. + /// + /// Returns the [`SymbolicAngle`] handles in order, so `angles[k]` is `α_k`. Allocating all of + /// a circuit's angles up front and then referring to them by index keeps the correspondence + /// between two circuits explicit and independent of how either circuit is otherwise written. + pub fn allocate_symbolic_angles(&mut self, count: usize) -> Vec { + (0..count).map(|_| self.allocate_symbolic_angle()).collect() + } + + /// All symbolic angles allocated on this simulation so far, in order (`angles[k]` is `α_k`). + #[getter] + #[must_use] + pub fn symbolic_angles(&self) -> Vec { + self.inner + .symbolic_angle_indicator() + .iter() + .enumerate() + .filter(|(_, &is_angle)| is_angle) + .enumerate() + .map(|(index, (outcome, _))| PySymbolicAngle { outcome, index }) + .collect() + } + + /// Apply a symbolic Pauli exponent `e^{iα P}` parameterised by `angle`. + /// + /// `angle` must be a [`SymbolicAngle`] obtained from [`allocate_symbolic_angle`] or + /// [`allocate_symbolic_angles`]. This is the high-level way to add a free-angle exponent + /// `e^{iα P}` for an arbitrary Pauli `P`. The same `angle` may parameterise several exponents + /// to model a shared `α`, and angles with matching `index` in two circuits are what make those + /// circuits' exponents correspond when their phased actions are compared. + #[allow(clippy::needless_pass_by_value)] + pub fn apply_symbolic_pauli_exp(&mut self, observable: &PySparsePauli, angle: &PySymbolicAngle) { + self.inner.symbolic_pauli_exp(&observable.inner, angle.outcome); + } + + #[allow(clippy::needless_pass_by_value)] + /// # Errors + /// + /// Returns a `ValueError` if the non-output system qubits remain entangled. + pub fn phased_action( + &self, + input_qubits: Vec, + output_qubits: Vec, + ) -> PyResult { + phased_action_from_simulation(&self.inner, &input_qubits, &output_qubits) + .map(|action| PyPhasedCircuitAction { inner: action }) + .map_err(|error| PyValueError::new_err(format!("{error:?}"))) + } +}); + +#[derive(derive_more::From)] +#[must_use] +#[pyclass(name = "PhasedCircuitAction", module = "paulimer")] +pub struct PyPhasedCircuitAction { + inner: PhasedCircuitAction, +} + +#[pymethods] +impl PyPhasedCircuitAction { + #[getter] + #[must_use] + pub fn choi_state_stabilizers(&self) -> Vec { + self.inner + .choi_state_stabilizers() + .iter() + .map(|pauli| PySparsePauli { inner: pauli.clone() }) + .collect() + } + + #[must_use] + pub fn is_equivalent(&self, other: &PyPhasedCircuitAction) -> bool { + self.inner.is_equivalent(&other.inner).is_ok() + } + + #[must_use] + pub fn is_equivalent_up_to_signs(&self, other: &PyPhasedCircuitAction) -> bool { + self.inner.is_equivalent_up_to_signs(&other.inner).is_ok() + } +} diff --git a/paulimer/bindings/python/tests/clifford_test.py b/paulimer/bindings/python/tests/clifford_test.py index 8cfb5219..411ced6e 100644 --- a/paulimer/bindings/python/tests/clifford_test.py +++ b/paulimer/bindings/python/tests/clifford_test.py @@ -320,6 +320,42 @@ def test_left_mul_pauli_exp_with_sparse_pauli(): assert clifford_dense.image_z(qubit) == clifford_sparse.image_z(qubit) +def _rebuild_from_pauli_exponents(exponents, num_qubits): + rebuilt = CliffordUnitary.identity(num_qubits) + for pauli in exponents: + rebuilt.left_mul_pauli_exp(pauli) + return rebuilt + + +def test_to_pauli_exponents_identity_is_empty(): + assert CliffordUnitary.identity(3).to_pauli_exponents() == [] + + +def test_to_pauli_exponents_roundtrip(): + clifford = CliffordUnitary.identity(3) + clifford.left_mul(UnitaryOpcode.Hadamard, [0]) + clifford.left_mul(UnitaryOpcode.ControlledX, [0, 1]) + clifford.left_mul(UnitaryOpcode.SqrtZ, [2]) + clifford.left_mul(UnitaryOpcode.ControlledZ, [1, 2]) + + exponents = clifford.to_pauli_exponents() + + assert all(isinstance(pauli, SparsePauli) for pauli in exponents) + assert _rebuild_from_pauli_exponents(exponents, 3) == clifford + + +def test_to_pauli_exponents_roundtrip_named_gates(): + for name, support, num_qubits in [ + ("Hadamard", [0], 1), + ("SqrtX", [0], 1), + ("ControlledX", [0, 1], 2), + ("ControlledZ", [0, 1], 2), + ]: + clifford = CliffordUnitary.from_name(name, support, num_qubits) + rebuilt = _rebuild_from_pauli_exponents(clifford.to_pauli_exponents(), num_qubits) + assert rebuilt == clifford + + def test_left_mul_controlled_pauli_with_dense_paulis(): clifford = CliffordUnitary.identity(2) control = DensePauli("ZI") diff --git a/paulimer/bindings/python/tests/simulation_test.py b/paulimer/bindings/python/tests/simulation_test.py index cacda7ad..e720771f 100644 --- a/paulimer/bindings/python/tests/simulation_test.py +++ b/paulimer/bindings/python/tests/simulation_test.py @@ -1,5 +1,7 @@ import pytest +import random from binar import BitMatrix, BitVector +from hypothesis import given, strategies as st from paulimer import ( CliffordUnitary, SparsePauli, @@ -7,12 +9,24 @@ OutcomeCompleteSimulation, OutcomeFreeSimulation, OutcomeSpecificSimulation, + PhasedCircuitAction, + PhasedOutcomeCompleteSimulation, ) SIMULATION_CLASSES = [ OutcomeCompleteSimulation, OutcomeFreeSimulation, OutcomeSpecificSimulation, + PhasedOutcomeCompleteSimulation, +] + +# PhasedOutcomeCompleteSimulation tracks the exact global phase, which a phaseless CliffordUnitary +# does not determine, so apply_clifford raises NotImplementedError on it; exclude it from the +# apply_clifford tests below while keeping it in the shared list for everything else. +CLIFFORD_CAPABLE_SIMULATION_CLASSES = [ + OutcomeCompleteSimulation, + OutcomeFreeSimulation, + OutcomeSpecificSimulation, ] @@ -130,13 +144,13 @@ def test_apply_permutation_full(self, sim_class): sim = sim_class(3) sim.apply_permutation([2, 0, 1]) - @pytest.mark.parametrize("sim_class", SIMULATION_CLASSES) + @pytest.mark.parametrize("sim_class", CLIFFORD_CAPABLE_SIMULATION_CLASSES) def test_apply_clifford_with_support(self, sim_class): sim = sim_class(3) hadamard = CliffordUnitary.from_name("Hadamard", [0], 1) sim.apply_clifford(hadamard, supported_by=[1]) - @pytest.mark.parametrize("sim_class", SIMULATION_CLASSES) + @pytest.mark.parametrize("sim_class", CLIFFORD_CAPABLE_SIMULATION_CLASSES) def test_apply_clifford_full(self, sim_class): sim = sim_class(2) cnot = CliffordUnitary.from_name("ControlledX", [0, 1], 2) @@ -258,3 +272,223 @@ def test_with_zero_outcomes(self): def test_new_with_seeded_random_outcomes(self): sim = OutcomeSpecificSimulation.new_with_seeded_random_outcomes(3, seed=42) assert sim.qubit_count == 3 + + +class TestPhasedOutcomeCompleteSimulationSpecific: + + def test_default_construction(self): + sim = PhasedOutcomeCompleteSimulation() + assert isinstance(sim.qubit_count, int) + + def test_construction_with_qubit_count(self): + sim = PhasedOutcomeCompleteSimulation(4) + assert sim.qubit_count == 4 + assert sim.qubit_capacity >= 4 + + def test_with_capacity(self): + sim = PhasedOutcomeCompleteSimulation.with_capacity(3, 10, 5) + assert sim.outcome_capacity >= 10 + assert sim.random_outcome_capacity >= 5 + + def test_clifford_returns_clifford_unitary(self): + sim = PhasedOutcomeCompleteSimulation(2) + assert isinstance(sim.clifford, CliffordUnitary) + + def test_phase_matrices_return_expected_types(self): + sim = PhasedOutcomeCompleteSimulation(2) + assert isinstance(sim.sign_matrix, BitMatrix) + assert isinstance(sim.quadratic_phase_matrix, BitMatrix) + assert isinstance(sim.outcome_matrix, BitMatrix) + assert isinstance(sim.outcome_shift, BitVector) + assert isinstance(sim.linear_i_phase, BitVector) + assert isinstance(sim.linear_sign_phase, BitVector) + + def test_operations_and_measurement(self): + sim = PhasedOutcomeCompleteSimulation(2) + sim.apply_unitary(UnitaryOpcode.Hadamard, [0]) + sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + outcome = sim.measure(SparsePauli("X_0")) + assert isinstance(outcome, int) + assert sim.random_outcome_count == 1 + + def test_output_phase_exponent(self): + sim = PhasedOutcomeCompleteSimulation(1) + sim.measure(SparsePauli("X_0")) + exponent = sim.output_phase_exponent([True]) + assert isinstance(exponent, int) + assert 0 <= exponent < 8 + # The trivial assignment never contributes a phase. + assert sim.output_phase_exponent([False]) == 0 + +def _choi_action(build_gadget, n=1): + """Phased Choi action of a symbolic-rotation gadget on ``n`` system qubits. + + Bell-pairs every system qubit ``q`` in ``0..n`` with its reference ``q + n``, allocates + one symbolic-angle random bit, applies the gadget to the system qubits, and returns the + resulting :class:`PhasedCircuitAction`. + """ + sim = PhasedOutcomeCompleteSimulation(2 * n) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.PrepareBell, [q, q + n]) + angle = sim.allocate_symbolic_angle() + build_gadget(sim, angle) + return sim.phased_action(list(range(n)), list(range(n))) + + +class TestPhasedCircuitAction: + + def test_phased_action_returns_action(self): + action = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_0"), a)) + assert isinstance(action, PhasedCircuitAction) + + def test_choi_state_stabilizers_are_sparse_paulis(self): + action = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_0"), a)) + stabilizers = action.choi_state_stabilizers + assert len(stabilizers) == 2 + assert all(isinstance(stabilizer, SparsePauli) for stabilizer in stabilizers) + + def test_entangling_rotation_equivalence(self): + def zz_direct(sim, a): + sim.apply_symbolic_pauli_exp(SparsePauli("Z_0 Z_1"), a) + + def zz_via_cnot(sim, a): + sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + sim.apply_symbolic_pauli_exp(SparsePauli("Z_1"), a) + sim.apply_unitary(UnitaryOpcode.ControlledX, [0, 1]) + + direct = _choi_action(zz_direct, n=2) + via_cnot = _choi_action(zz_via_cnot, n=2) + assert direct.is_equivalent(via_cnot) + assert via_cnot.is_equivalent(direct) + + def test_dropping_conjugation_breaks_equivalence(self): + direct = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_0 Z_1"), a), n=2) + bare = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_1"), a), n=2) + assert not direct.is_equivalent(bare) + + def test_opposite_signs_distinguished_only_by_phase(self): + positive = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_0"), a)) + negative = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("-Z_0"), a)) + assert positive.is_equivalent_up_to_signs(negative) + assert not positive.is_equivalent(negative) + + def test_action_is_self_equivalent(self): + action = _choi_action(lambda sim, a: sim.apply_symbolic_pauli_exp(SparsePauli("Z_0 Z_1"), a), n=2) + assert action.is_equivalent(action) + + +def _direct_z_action(n, angle_supports): + """Phased action of symbolic Z-rotations applied directly to ``n`` system qubits. + + Layout: system qubits ``0..n``, reference qubits ``n..2n``. Each entry of + ``angle_supports`` is a list of system-qubit indices naming one symbolic rotation + ``e^{i alpha Z...}`` about the tensor product of ``Z`` on those qubits. + """ + sim = PhasedOutcomeCompleteSimulation(2 * n) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.PrepareBell, [q, n + q]) + for support in angle_supports: + pauli = SparsePauli(" ".join(f"Z_{q}" for q in support)) + angle = sim.allocate_symbolic_angle() + sim.apply_symbolic_pauli_exp(pauli, angle) + return sim.phased_action(list(range(n)), list(range(n))) + + +def _z_ejection_action(n, angle_supports): + """Phased action of the same symbolic Z-rotations executed remotely via ejection. + + Layout: system qubits ``0..n``, reference qubits ``n..2n``, ancillas ``2n..3n``. The + system qubits drive transversal CNOTs onto the ancillas, the symbolic rotations act on + the ancillas, each ancilla is measured in the X basis (a *true* random bit), and a + ``-`` outcome triggers a conditional Z correction on the system qubit. Mixing the virtual + angle bits with the true measurement bits, the action must equal :func:`_direct_z_action`. + """ + sim = PhasedOutcomeCompleteSimulation(3 * n) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.PrepareBell, [q, n + q]) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.ControlledX, [q, 2 * n + q]) + for support in angle_supports: + pauli = SparsePauli(" ".join(f"Z_{2 * n + q}" for q in support)) + angle = sim.allocate_symbolic_angle() + sim.apply_symbolic_pauli_exp(pauli, angle) + for q in range(n): + outcome = sim.measure(SparsePauli(f"X_{2 * n + q}")) + sim.apply_conditional_pauli(SparsePauli(f"Z_{q}"), [outcome]) + return sim.phased_action(list(range(n)), list(range(n))) + + +class TestPhasedEjection: + """Measurement-based "ejection" of a Z-diagonal channel, mixing virtual and true bits.""" + + @pytest.mark.parametrize( + "n, angle_supports", + [ + (1, [[0]]), + (2, [[0]]), + (2, [[0, 1]]), + (2, [[0], [1], [0, 1]]), + ], + ) + def test_ejection_matches_direct(self, n, angle_supports): + direct = _direct_z_action(n, angle_supports) + ejection = _z_ejection_action(n, angle_supports) + assert direct.is_equivalent(ejection) + assert ejection.is_equivalent(direct) + +def _distinct_z_supports(n): + """All non-trivial Z-product supports on ``n`` qubits (distinct, independent).""" + return [[bit for bit in range(n) if mask & (1 << bit)] for mask in range(1, 1 << n)] + + +def _signed_direct_z_action(n, supports, signs): + """:func:`_direct_z_action` where the ``k``-th rotation is negated iff ``signs[k]``.""" + sim = PhasedOutcomeCompleteSimulation(2 * n) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.PrepareBell, [q, n + q]) + for support, negate in zip(supports, signs): + body = " ".join(f"Z_{q}" for q in support) + pauli = SparsePauli(("-" if negate else "") + body) + sim.apply_symbolic_pauli_exp(pauli, sim.allocate_symbolic_angle()) + return sim.phased_action(list(range(n)), list(range(n))) + + +def _permuted_direct_z_action(n, supports, perm): + """:func:`_direct_z_action` where allocated angle ``k`` drives ``supports[perm[k]]``.""" + sim = PhasedOutcomeCompleteSimulation(2 * n) + for q in range(n): + sim.apply_unitary(UnitaryOpcode.PrepareBell, [q, n + q]) + for target in perm: + pauli = SparsePauli(" ".join(f"Z_{q}" for q in supports[target])) + sim.apply_symbolic_pauli_exp(pauli, sim.allocate_symbolic_angle()) + return sim.phased_action(list(range(n)), list(range(n))) + + +class TestPhasedNegativeChecks: + """Sign flips and angle permutations must break an otherwise exact phased equivalence.""" + + @given( + n=st.integers(min_value=2, max_value=3), + seed=st.integers(min_value=0, max_value=2**32 - 1), + ) + def test_sign_flip_is_pure_relative_phase(self, n, seed): + supports = _distinct_z_supports(n) + rng = random.Random(seed) + signs = [rng.random() < 0.5 for _ in supports] + baseline = _signed_direct_z_action(n, supports, [False] * len(supports)) + flipped = _signed_direct_z_action(n, supports, signs) + assert baseline.is_equivalent_up_to_signs(flipped) + assert baseline.is_equivalent(flipped) == (not any(signs)) + + @given( + n=st.integers(min_value=2, max_value=3), + seed=st.integers(min_value=0, max_value=2**32 - 1), + ) + def test_permuted_angles_match_iff_identity(self, n, seed): + supports = _distinct_z_supports(n) + identity = list(range(len(supports))) + perm = identity[:] + random.Random(seed).shuffle(perm) + base = _permuted_direct_z_action(n, supports, identity) + permuted = _permuted_direct_z_action(n, supports, perm) + assert base.is_equivalent(permuted) == (perm == identity) diff --git a/paulimer/src/clifford.rs b/paulimer/src/clifford.rs index 5fa75826..bdac6257 100644 --- a/paulimer/src/clifford.rs +++ b/paulimer/src/clifford.rs @@ -299,6 +299,8 @@ pub struct CliffordModPauliBatch Vec { + let mut recorded = Reduction::new(); + let mut working = clifford.clone(); + for pivot in 0..working.num_qubits() { + recorded.clear_x_image(&mut working, pivot); + recorded.clear_z_image(&mut working, pivot); + } + debug_assert!(working.is_identity(), "Clifford reduction did not reach the identity"); + recorded.into_decomposition() +} + +/// Accumulates the `π/4` exponents applied while reducing a Clifford to the identity. +struct Reduction { + applied: Vec, +} + +impl Reduction { + fn new() -> Self { + Reduction { applied: Vec::new() } + } + + /// Left-multiplies `working` by `exp(iπ/4·pauli)` and records the factor. + fn exp(&mut self, working: &mut CliffordUnitary, pauli: SparsePauli) { + working.left_mul_pauli_exp(&pauli); + self.applied.push(pauli); + } + + fn single_x(qubit: usize, qubit_count: usize) -> SparsePauli { + SparsePauli::x(qubit, qubit_count) + } + + fn single_z(qubit: usize, qubit_count: usize) -> SparsePauli { + SparsePauli::z(qubit, qubit_count) + } + + /// `Z_control · X_target`, the generator of a controlled-`X`. + fn control_x(control: usize, target: usize, qubit_count: usize) -> SparsePauli { + let mut pauli = SparsePauli::z(control, qubit_count); + pauli.mul_assign_left_x(target); + pauli + } + + /// `Z_a · Z_b`, the generator of a controlled-`Z`. + fn control_z(first: usize, second: usize, qubit_count: usize) -> SparsePauli { + let mut pauli = SparsePauli::z(first, qubit_count); + pauli.mul_assign_left_z(second); + pauli + } + + fn hadamard(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, Self::single_x(qubit, count)); + self.exp(working, Self::single_z(qubit, count)); + self.exp(working, Self::single_x(qubit, count)); + } + + fn root_z(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, Self::single_z(qubit, count)); + } + + fn root_z_inverse(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, negated(Self::single_z(qubit, count))); + } + + fn root_x(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, Self::single_x(qubit, count)); + } + + fn controlled_x(&mut self, working: &mut CliffordUnitary, control: usize, target: usize) { + let count = working.num_qubits(); + self.exp(working, Self::control_x(control, target, count)); + self.exp(working, negated(Self::single_z(control, count))); + self.exp(working, negated(Self::single_x(target, count))); + } + + fn controlled_z(&mut self, working: &mut CliffordUnitary, first: usize, second: usize) { + let count = working.num_qubits(); + self.exp(working, Self::control_z(first, second, count)); + self.exp(working, negated(Self::single_z(first, count))); + self.exp(working, negated(Self::single_z(second, count))); + } + + /// Conjugation by the Pauli `Z_qubit`, flipping the sign of an `X`-type image on `qubit`. + fn pauli_z(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, Self::single_z(qubit, count)); + self.exp(working, Self::single_z(qubit, count)); + } + + /// Conjugation by the Pauli `X_qubit`, flipping the sign of a `Z`-type image on `qubit`. + fn pauli_x(&mut self, working: &mut CliffordUnitary, qubit: usize) { + let count = working.num_qubits(); + self.exp(working, Self::single_x(qubit, count)); + self.exp(working, Self::single_x(qubit, count)); + } + + /// Turns the image of `X_pivot` into `+X_pivot` using gates supported on qubits `≥ pivot`. + fn clear_x_image(&mut self, working: &mut CliffordUnitary, pivot: usize) { + let count = working.num_qubits(); + let image = working.image_x(pivot); + if !(pivot..count).any(|qubit| x_bit(&image, qubit)) { + let qubit = (pivot..count) + .find(|&qubit| z_bit(&image, qubit)) + .expect("a non-identity image has an X or Z component"); + self.hadamard(working, qubit); + } + let image = working.image_x(pivot); + if !x_bit(&image, pivot) { + let qubit = (pivot..count) + .find(|&qubit| qubit != pivot && x_bit(&image, qubit)) + .expect("the image has an X component to move onto the pivot"); + self.controlled_x(working, qubit, pivot); + } + let image = working.image_x(pivot); + for qubit in (pivot + 1)..count { + if x_bit(&image, qubit) { + self.controlled_x(working, pivot, qubit); + } + } + let image = working.image_x(pivot); + if z_bit(&image, pivot) { + self.root_z_inverse(working, pivot); + } + let image = working.image_x(pivot); + for qubit in (pivot + 1)..count { + if z_bit(&image, qubit) { + self.controlled_z(working, pivot, qubit); + } + } + if working.image_x(pivot).xz_phase_exponent() != 0 { + self.pauli_z(working, pivot); + } + } + + /// Turns the image of `Z_pivot` into `+Z_pivot`, assuming the image of `X_pivot` is already + /// `+X_pivot`; every gate used fixes `X_pivot`. + fn clear_z_image(&mut self, working: &mut CliffordUnitary, pivot: usize) { + let count = working.num_qubits(); + if x_bit(&working.image_z(pivot), pivot) { + self.root_x(working, pivot); + } + for qubit in (pivot + 1)..count { + let image = working.image_z(pivot); + if x_bit(&image, qubit) && z_bit(&image, qubit) { + self.root_z(working, qubit); + } + if x_bit(&working.image_z(pivot), qubit) { + self.hadamard(working, qubit); + } + if z_bit(&working.image_z(pivot), qubit) { + self.controlled_x(working, qubit, pivot); + } + } + if working.image_z(pivot).xz_phase_exponent() != 0 { + self.pauli_x(working, pivot); + } + } + + /// The exponents that rebuild the original Clifford from the identity. + fn into_decomposition(self) -> Vec { + self.applied.into_iter().rev().map(negated).collect() + } +} + +/// `exp(iπ/4·P)⁻¹ = exp(iπ/4·(−P))`, with `−P` encoded as a phase-exponent shift of two. +fn negated(mut pauli: SparsePauli) -> SparsePauli { + pauli.add_assign_phase_exp(2); + pauli +} + +fn x_bit(pauli: &P, qubit: usize) -> bool { + pauli.x_bits().index(qubit) +} + +fn z_bit(pauli: &P, qubit: usize) -> bool { + pauli.z_bits().index(qubit) +} diff --git a/paulimer/src/clifford/phased_clifford.rs b/paulimer/src/clifford/phased_clifford.rs new file mode 100644 index 00000000..2c499215 --- /dev/null +++ b/paulimer/src/clifford/phased_clifford.rs @@ -0,0 +1,479 @@ +//! Phase-tracking Clifford unitaries. +//! +//! A [`CliffordUnitary`] represents a Clifford operator only up to a global phase: its symplectic +//! matrix together with the signs of the Pauli images fixes the operator on the Pauli group, but +//! not the overall `ζ₈` factor of the unitary. Many stabilizer algorithms do not need this factor, +//! but the *phased* outcome-complete simulation of arXiv:2603.24717 does, because it tracks the +//! exact amplitudes (including global phase) of the simulated state. +//! +//! [`PhasedCliffordUnitary`] augments a [`CliffordUnitary`] with an exact global-phase tracker. +//! The key observation is that everywhere a phase is needed it is the phase of the *encoder state* +//! `C|0…0⟩`: for any bit string `a`, `C|a⟩ = C X^a C† · C|0…0⟩`, and the sign of the Pauli image +//! `C X^a C†` is already tracked by [`CliffordUnitary`]. It is therefore enough to maintain one +//! exactly-known amplitude of `C|0…0⟩`, namely its value at a fixed *reference* basis string. +//! +//! Concretely the tracker stores a basis string `r` with `⟨r|C|0…0⟩ ≠ 0` and the `ζ₈` exponent of +//! that amplitude. Every amplitude of the stabilizer state `C|0…0⟩` has the same magnitude, so the +//! magnitude is recovered from the rank of the stabilizer tableau and only the phase needs to be +//! propagated. Each elementary left-multiplication updates the underlying [`CliffordUnitary`] with +//! the existing tableau code and updates the reference amplitude in `O(n²)` time. + +use super::{Clifford, CliffordMutable, CliffordUnitary}; +use crate::UnitaryOp; +use crate::pauli::{Pauli, PauliBinaryOps}; +use binar::matrix::AlignedBitMatrix; +use binar::vec::AlignedBitVec; +use binar::{BitMatrix, BitVec, Bitwise, BitwiseMut, EchelonForm}; + +fn normalize_exponent(value: i64) -> u8 { + u8::try_from(value.rem_euclid(8)).expect("a value reduced modulo 8 is in 0..8") +} + +/// Sums two unit `ζ₈` powers `ζ₈^x + ζ₈^y` exactly, returning the `ζ₈` exponent of the result. +/// +/// In this simulator every such sum that arises is the amplitude of a stabilizer state, hence a +/// non-negative real multiple of a single `ζ₈` power. Two unit `ζ₈` powers can therefore only be +/// separated by `d = (x − y) mod 8 ∈ {0, 2, 6}` (parallel or orthogonal), or by `d = 4`, in which +/// case they cancel exactly. This makes the sum computable with pure integer arithmetic, with no +/// floating-point amplitudes. Returns `None` when the two terms cancel. +fn add_zeta8_powers(x: i64, y: i64) -> Option { + match (x - y).rem_euclid(8) { + 0 => Some(normalize_exponent(x)), + 2 => Some(normalize_exponent(y + 1)), + 6 => Some(normalize_exponent(y + 7)), + 4 => None, + separation => unreachable!("ζ₈^{x} + ζ₈^{y} is not a ζ₈ power (separation {separation})"), + } +} + +/// Combines the at-most-two unit `ζ₈` powers contributing to a single amplitude, returning the +/// `ζ₈` exponent of their sum, or `None` when the amplitude vanishes. See [`add_zeta8_powers`]. +fn combine_zeta8_powers(terms: &[i64]) -> Option { + match *terms { + [] => None, + [only] => Some(normalize_exponent(only)), + [x, y] => add_zeta8_powers(x, y), + _ => unreachable!("an encoder amplitude is a sum of at most two ζ₈ powers"), + } +} + +/// A Clifford unitary that additionally tracks the exact global phase of its encoder state. +/// +/// The operator agrees with [`Self::clifford`] on the Pauli group (same symplectic action and same +/// image signs) and in addition fixes the overall `ζ₈` phase, so that the amplitudes of `C|0…0⟩` +/// are determined exactly rather than only up to a global factor. +/// +/// # Examples +/// +/// ``` +/// use paulimer::clifford::PhasedCliffordUnitary; +/// +/// let mut phased = PhasedCliffordUnitary::identity(1); +/// phased.left_mul_hadamard(0); +/// // |0⟩ -> (|0⟩ + |1⟩)/√2: the amplitude at |0⟩ is positive real (ζ₈ exponent 0). +/// assert_eq!(phased.state_amplitude_phase_exponent_usize(0), Some(0)); +/// ``` +#[must_use] +#[derive(Clone)] +pub struct PhasedCliffordUnitary { + clifford: CliffordUnitary, + reference_string: AlignedBitVec, + reference_phase_exponent: u8, +} + +impl PhasedCliffordUnitary { + /// Returns the identity operator on `num_qubits` qubits, with encoder state `|0…0⟩`. + pub fn identity(num_qubits: usize) -> Self { + Self { + clifford: CliffordUnitary::identity(num_qubits), + reference_string: AlignedBitVec::zeros(num_qubits), + reference_phase_exponent: 0, + } + } + + /// Returns the number of qubits the operator acts on. + #[must_use] + pub fn num_qubits(&self) -> usize { + self.clifford.num_qubits() + } + + /// Returns the underlying phaseless [`CliffordUnitary`]. + pub fn clifford(&self) -> &CliffordUnitary { + &self.clifford + } + + /// Consumes the operator and returns the underlying phaseless [`CliffordUnitary`]. + pub fn into_clifford(self) -> CliffordUnitary { + self.clifford + } + + /// Left-multiplies by the global scalar `ζ₈^exponent = e^{i π exponent / 4}`. + /// + /// This leaves the underlying [`CliffordUnitary`] (and hence the symplectic action and image + /// signs) unchanged and only advances the tracked global phase, since multiplying the operator + /// by a scalar scales every amplitude of the encoder state by the same factor. + pub fn left_mul_global_phase(&mut self, exponent: u8) { + self.reference_phase_exponent = + normalize_exponent(i64::from(self.reference_phase_exponent) + i64::from(exponent)); + } + + /// Returns the `ζ₈` exponent `e` such that `⟨basis|C|0…0⟩ = ζ₈^e · 2^{-k/2}` for some rank `k`, + /// or `None` when that amplitude vanishes. + #[must_use] + pub fn state_amplitude_phase_exponent(&self, basis: &AlignedBitVec) -> Option { + let relative = self.relative_phase(basis)?; + Some(normalize_exponent(i64::from(self.reference_phase_exponent) + relative)) + } + + /// Convenience wrapper around [`Self::state_amplitude_phase_exponent`] taking the basis string + /// as an integer whose qubit `q` bit is `(value >> q) & 1`. + #[must_use] + pub fn state_amplitude_phase_exponent_usize(&self, value: usize) -> Option { + let mut basis = AlignedBitVec::zeros(self.num_qubits()); + for qubit in 0..self.num_qubits() { + basis.assign_index(qubit, (value >> qubit) & 1 == 1); + } + self.state_amplitude_phase_exponent(&basis) + } + + fn x_parts_matrix(&self) -> BitMatrix { + let num_qubits = self.num_qubits(); + let mut matrix = AlignedBitMatrix::zeros(num_qubits, num_qubits); + for generator in 0..num_qubits { + let image = self.clifford.image_z(generator); + for qubit in image.x_bits().support() { + matrix.row_mut(generator).assign_index(qubit, true); + } + } + BitMatrix::from_aligned(matrix) + } + + fn relative_phase(&self, target: &AlignedBitVec) -> Option { + let num_qubits = self.num_qubits(); + let difference: BitVec = (0..num_qubits) + .map(|qubit| target.index(qubit) ^ self.reference_string.index(qubit)) + .collect(); + let echelon = EchelonForm::new(self.x_parts_matrix()); + let combination = echelon.transpose_solve(&difference.as_view())?; + let mut product = self.clifford.image_z(0); + let mut started = false; + for generator in combination.support() { + let image = self.clifford.image_z(generator); + if started { + product.mul_assign_right(&image); + } else { + product = image; + started = true; + } + } + if !started { + return Some(0); + } + let phase_exponent = i64::from(product.xz_phase_exponent()); + let mut sign_parity = false; + for qubit in product.z_bits().support() { + if self.reference_string.index(qubit) { + sign_parity = !sign_parity; + } + } + let relative = (2 * phase_exponent + if sign_parity { 4 } else { 0 }).rem_euclid(8); + Some(relative) + } + + fn apply_one_qubit( + &mut self, + qubit: usize, + amplitudes: [[Option; 2]; 2], + symplectic: impl FnOnce(&mut CliffordUnitary), + ) { + for output_bit in [self.reference_string.index(qubit), !self.reference_string.index(qubit)] { + let mut candidate = self.reference_string.clone(); + candidate.assign_index(qubit, output_bit); + let mut terms = [0i64; 2]; + let mut count = 0usize; + for input_bit in [false, true] { + let Some(entry) = amplitudes[usize::from(output_bit)][usize::from(input_bit)] else { + continue; + }; + let mut source = candidate.clone(); + source.assign_index(qubit, input_bit); + if let Some(relative) = self.relative_phase(&source) { + terms[count] = entry + relative; + count += 1; + } + } + if let Some(increment) = combine_zeta8_powers(&terms[..count]) { + self.reference_phase_exponent = + normalize_exponent(i64::from(self.reference_phase_exponent) + i64::from(increment)); + self.reference_string = candidate; + symplectic(&mut self.clifford); + return; + } + } + unreachable!("a unitary maps a nonzero state to a nonzero state"); + } + + fn apply_two_qubit( + &mut self, + qubit_a: usize, + qubit_b: usize, + inverse: impl Fn(bool, bool) -> (bool, bool, i64), + symplectic: impl FnOnce(&mut CliffordUnitary), + ) { + for output_a in [ + self.reference_string.index(qubit_a), + !self.reference_string.index(qubit_a), + ] { + for output_b in [ + self.reference_string.index(qubit_b), + !self.reference_string.index(qubit_b), + ] { + let mut candidate = self.reference_string.clone(); + candidate.assign_index(qubit_a, output_a); + candidate.assign_index(qubit_b, output_b); + let (input_a, input_b, entry) = inverse(output_a, output_b); + let mut source = candidate.clone(); + source.assign_index(qubit_a, input_a); + source.assign_index(qubit_b, input_b); + if let Some(relative) = self.relative_phase(&source) { + self.reference_phase_exponent = + normalize_exponent(i64::from(self.reference_phase_exponent) + entry + relative); + self.reference_string = candidate; + symplectic(&mut self.clifford); + return; + } + } + } + unreachable!("a unitary maps a nonzero state to a nonzero state"); + } + + /// Left-multiplies by a Hadamard gate on `qubit`. + pub fn left_mul_hadamard(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(0), Some(0)], [Some(0), Some(4)]], |clifford| { + clifford.left_mul_hadamard(qubit); + }); + } + + /// Left-multiplies by a Pauli `X` gate on `qubit`. + pub fn left_mul_x(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[None, Some(0)], [Some(0), None]], |clifford| { + clifford.left_mul_x(qubit); + }); + } + + /// Left-multiplies by a Pauli `Y` gate on `qubit`. + pub fn left_mul_y(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[None, Some(6)], [Some(2), None]], |clifford| { + clifford.left_mul_y(qubit); + }); + } + + /// Left-multiplies by a Pauli `Z` gate on `qubit`. + pub fn left_mul_z(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(0), None], [None, Some(4)]], |clifford| { + clifford.left_mul_z(qubit); + }); + } + + /// Left-multiplies by `√Z` (the phase gate `S = diag(1, i)`) on `qubit`. + pub fn left_mul_root_z(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(0), None], [None, Some(2)]], |clifford| { + clifford.left_mul_root_z(qubit); + }); + } + + /// Left-multiplies by `√Z†` (`S† = diag(1, -i)`) on `qubit`. + pub fn left_mul_root_z_inverse(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(0), None], [None, Some(6)]], |clifford| { + clifford.left_mul_root_z_inverse(qubit); + }); + } + + /// Left-multiplies by `√X` on `qubit`. + pub fn left_mul_root_x(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(1), Some(7)], [Some(7), Some(1)]], |clifford| { + clifford.left_mul_root_x(qubit); + }); + } + + /// Left-multiplies by `√X†` on `qubit`. + pub fn left_mul_root_x_inverse(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(7), Some(1)], [Some(1), Some(7)]], |clifford| { + clifford.left_mul_root_x_inverse(qubit); + }); + } + + /// Left-multiplies by `√Y` on `qubit`. + pub fn left_mul_root_y(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(1), Some(5)], [Some(1), Some(1)]], |clifford| { + clifford.left_mul_root_y(qubit); + }); + } + + /// Left-multiplies by `√Y†` on `qubit`. + pub fn left_mul_root_y_inverse(&mut self, qubit: usize) { + self.apply_one_qubit(qubit, [[Some(7), Some(7)], [Some(3), Some(7)]], |clifford| { + clifford.left_mul_root_y_inverse(qubit); + }); + } + + /// Left-multiplies by a controlled-`X` gate with the given control and target qubits. + pub fn left_mul_cx(&mut self, control: usize, target: usize) { + self.apply_two_qubit( + control, + target, + |control_bit, target_bit| (control_bit, control_bit ^ target_bit, 0), + |clifford| { + clifford.left_mul_cx(control, target); + }, + ); + } + + /// Left-multiplies by a controlled-`Z` gate on the two given qubits. + pub fn left_mul_cz(&mut self, qubit_a: usize, qubit_b: usize) { + self.apply_two_qubit( + qubit_a, + qubit_b, + |bit_a, bit_b| (bit_a, bit_b, if bit_a && bit_b { 4 } else { 0 }), + |clifford| { + clifford.left_mul_cz(qubit_a, qubit_b); + }, + ); + } + + /// Left-multiplies by a swap of the two given qubits. + pub fn left_mul_swap(&mut self, qubit_a: usize, qubit_b: usize) { + self.apply_two_qubit( + qubit_a, + qubit_b, + |bit_a, bit_b| (bit_b, bit_a, 0), + |clifford| { + clifford.left_mul_swap(qubit_a, qubit_b); + }, + ); + } + + /// Left-multiplies by the named elementary [`UnitaryOp`] on `support`. + pub fn left_mul(&mut self, unitary_op: UnitaryOp, support: &[usize]) { + use UnitaryOp::{ + ControlledX, ControlledZ, Hadamard, I, PrepareBell, SqrtX, SqrtXInv, SqrtY, SqrtYInv, SqrtZ, SqrtZInv, + Swap, X, Y, Z, + }; + match unitary_op { + I => {} + X => self.left_mul_x(support[0]), + Y => self.left_mul_y(support[0]), + Z => self.left_mul_z(support[0]), + SqrtX => self.left_mul_root_x(support[0]), + SqrtXInv => self.left_mul_root_x_inverse(support[0]), + SqrtY => self.left_mul_root_y(support[0]), + SqrtYInv => self.left_mul_root_y_inverse(support[0]), + SqrtZ => self.left_mul_root_z(support[0]), + SqrtZInv => self.left_mul_root_z_inverse(support[0]), + Hadamard => self.left_mul_hadamard(support[0]), + Swap => self.left_mul_swap(support[0], support[1]), + ControlledX => self.left_mul_cx(support[0], support[1]), + ControlledZ => self.left_mul_cz(support[0], support[1]), + PrepareBell => self.left_mul_prepare_bell(support[0], support[1]), + } + } + + /// Left-multiplies by the qubit permutation `permutation` acting on `support`. + /// + /// A permutation of the computational basis labels has trivial global phase, so only the + /// underlying tableau and the reference basis string are relabelled while the tracked phase is + /// left unchanged. The convention matches [`CliffordMutable::left_mul_permutation`]: the qubit + /// `support[i]` takes the role previously played by `support[permutation[i]]`. + pub fn left_mul_permutation(&mut self, permutation: &[usize], support: &[usize]) { + let previous: Vec = support + .iter() + .map(|&qubit| self.reference_string.index(qubit)) + .collect(); + self.clifford.left_mul_permutation(permutation, support); + for (index, &qubit) in support.iter().enumerate() { + self.reference_string.assign_index(qubit, previous[permutation[index]]); + } + } + + /// Left-multiplies by the Bell-state preparation Clifford on the two given qubits. + pub fn left_mul_prepare_bell(&mut self, qubit_a: usize, qubit_b: usize) { + self.left_mul_hadamard(qubit_a); + self.left_mul_cx(qubit_a, qubit_b); + } + + /// Left-multiplies by the Pauli operator `pauli` (including its sign). + pub fn left_mul_pauli>(&mut self, pauli: &PauliLike) { + let phase = pauli.xz_phase_exponent(); + for qubit in pauli.z_bits().support() { + self.left_mul_z(qubit); + } + for qubit in pauli.x_bits().support() { + self.left_mul_x(qubit); + } + if phase != 0 { + self.reference_phase_exponent = + normalize_exponent(i64::from(self.reference_phase_exponent) + 2 * i64::from(phase)); + } + } + + /// Left-multiplies by `exp(iπ/4 · pauli)`, the square root of `pauli` up to phase. + pub fn left_mul_pauli_exp>(&mut self, pauli: &PauliLike) { + if self.num_qubits() == 0 { + return; + } + let pauli_phase = i64::from(pauli.xz_phase_exponent()); + + let mut shifted = self.reference_string.clone(); + for qubit in pauli.x_bits().support() { + shifted.assign_index(qubit, !shifted.index(qubit)); + } + let candidates = [self.reference_string.clone(), shifted]; + + for candidate in candidates { + let mut terms = [0i64; 2]; + let mut count = 0usize; + if let Some(relative) = self.relative_phase(&candidate) { + terms[count] = relative; + count += 1; + } + let mut source = candidate.clone(); + let mut sign_parity = false; + for qubit in pauli.z_bits().support() { + let flipped = source.index(qubit) ^ pauli.x_bits().index(qubit); + if flipped { + sign_parity = !sign_parity; + } + } + for qubit in pauli.x_bits().support() { + source.assign_index(qubit, !source.index(qubit)); + } + if let Some(relative) = self.relative_phase(&source) { + let coefficient = 2 + 2 * pauli_phase + if sign_parity { 4 } else { 0 }; + terms[count] = relative + coefficient; + count += 1; + } + if let Some(increment) = combine_zeta8_powers(&terms[..count]) { + self.reference_phase_exponent = + normalize_exponent(i64::from(self.reference_phase_exponent) + i64::from(increment)); + self.reference_string = candidate; + self.clifford.left_mul_pauli_exp(pauli); + return; + } + } + unreachable!("a unitary maps a nonzero state to a nonzero state"); + } + + /// Grows or shrinks the operator to `new_num_qubits`, tensoring with identity on `|0⟩` ancillas. + pub fn resize(&mut self, new_num_qubits: usize) { + let old_num_qubits = self.num_qubits(); + self.clifford.resize(new_num_qubits); + if new_num_qubits == old_num_qubits { + return; + } + let mut reference = AlignedBitVec::zeros(new_num_qubits); + for qubit in 0..old_num_qubits.min(new_num_qubits) { + reference.assign_index(qubit, self.reference_string.index(qubit)); + } + self.reference_string = reference; + } +} diff --git a/paulimer/tests/clifford_test.rs b/paulimer/tests/clifford_test.rs index a45117d6..47d3e27f 100644 --- a/paulimer/tests/clifford_test.rs +++ b/paulimer/tests/clifford_test.rs @@ -7,8 +7,8 @@ use paulimer::StringNotation::{Ascii, Tex, Unicode}; use paulimer::clifford::generic_algos::{clifford_from_images, clifford_to_prepare_bell_states}; use paulimer::clifford::{ Clifford, CliffordMutable, CliffordStringParsingError, MutablePreImages, PreimageViews, XOrZ, - apply_qubit_clifford_by_axis, group_encoding_clifford_of, prepare_all_plus, prepare_all_zero, - random_clifford_via_operations_sampling, split_clifford_encoder_mod_pauli, split_phased_css, + apply_qubit_clifford_by_axis, clifford_to_pauli_exponents, group_encoding_clifford_of, prepare_all_plus, + prepare_all_zero, random_clifford_via_operations_sampling, split_clifford_encoder_mod_pauli, split_phased_css, split_qubit_cliffords_and_css, split_qubit_tensor_product_encoder, standard_restriction_with_sign_matrix, z_images_partition_transform, }; @@ -509,6 +509,43 @@ prop_compose! { } } +fn reconstruct_from_pauli_exponents(exponents: &[SparsePauli], dimension: usize) -> CliffordUnitary { + let mut rebuilt = CliffordUnitary::identity(dimension); + for exponent in exponents { + rebuilt.left_mul_pauli_exp(exponent); + } + rebuilt +} + +#[test] +fn clifford_to_pauli_exponents_identity_is_empty() { + for dimension in 0..4 { + let exponents = clifford_to_pauli_exponents(&CliffordUnitary::identity(dimension)); + assert!( + exponents.is_empty(), + "identity decomposes to no exponents (dimension {dimension})" + ); + } +} + +#[test] +fn clifford_to_pauli_exponents_examples_roundtrip() { + for clifford in clifford_examples::() { + let dimension = clifford.num_qubits(); + let exponents = clifford_to_pauli_exponents(&clifford); + assert_eq!(reconstruct_from_pauli_exponents(&exponents, dimension), clifford); + } +} + +proptest! { + #[test] + fn clifford_to_pauli_exponents_roundtrip(clifford in arbitrary_clifford(0..6)) { + let dimension = clifford.num_qubits(); + let exponents = clifford_to_pauli_exponents(&clifford); + prop_assert_eq!(reconstruct_from_pauli_exponents(&exponents, dimension), clifford); + } +} + prop_compose! { fn arbitrary_css_clifford(dimension_range: Range)(dimension in dimension_range) -> CliffordUnitary { let mut clifford: CliffordUnitary = random_css_clifford(dimension); diff --git a/paulimer/tests/phased_clifford_dense.rs b/paulimer/tests/phased_clifford_dense.rs new file mode 100644 index 00000000..5ee69ecd --- /dev/null +++ b/paulimer/tests/phased_clifford_dense.rs @@ -0,0 +1,102 @@ +use dense_oracle::{Dense, close, gate_matrix, pauli_arrays, statevector}; +use paulimer::clifford::PhasedCliffordUnitary; +use paulimer::pauli::Pauli; +use paulimer::{DensePauli, UnitaryOp}; +use proptest::collection::vec; +use proptest::prelude::*; + +#[derive(Clone, Debug)] +enum Gate { + Single { op: UnitaryOp, qubit: usize }, + Two { op: UnitaryOp, first: usize, second: usize }, + Pauli(DensePauli), + PauliExp(DensePauli), +} + +fn distinct_pair(qubit_count: usize) -> impl Strategy { + (0..qubit_count, 0..qubit_count - 1) + .prop_map(|(first, second)| (first, if second < first { second } else { second + 1 })) +} + +fn pauli_strategy(qubit_count: usize, hermitian: bool) -> impl Strategy { + (vec(any::(), qubit_count), vec(any::(), qubit_count), 0u8..4) + .prop_map(|(x_bits, z_bits, phase)| { + DensePauli::from_bits(x_bits.into_iter().collect(), z_bits.into_iter().collect(), phase) + }) + .prop_filter("non-identity, Hermitian when required", move |pauli| { + !pauli.is_identity() && (!hermitian || pauli.is_order_two()) + }) +} + +fn gate_strategy(qubit_count: usize) -> impl Strategy { + use UnitaryOp::{ + ControlledX, ControlledZ, Hadamard, PrepareBell, SqrtX, SqrtXInv, SqrtY, SqrtYInv, SqrtZ, SqrtZInv, Swap, X, Y, + Z, + }; + let single = ( + prop::sample::select(vec![ + Hadamard, X, Y, Z, SqrtZ, SqrtZInv, SqrtX, SqrtXInv, SqrtY, SqrtYInv, + ]), + 0..qubit_count, + ) + .prop_map(|(op, qubit)| Gate::Single { op, qubit }); + let two = ( + prop::sample::select(vec![ControlledX, ControlledZ, Swap, PrepareBell]), + distinct_pair(qubit_count), + ) + .prop_map(|(op, (first, second))| Gate::Two { op, first, second }); + prop_oneof![ + 10 => single, + 4 => two, + 1 => pauli_strategy(qubit_count, false).prop_map(Gate::Pauli), + 1 => pauli_strategy(qubit_count, true).prop_map(Gate::PauliExp), + ] +} + +fn apply(gate: &Gate, dense: &mut Dense, phased: &mut PhasedCliffordUnitary) { + use UnitaryOp::{ControlledX, ControlledZ, Hadamard, PrepareBell, Swap}; + let qubit_count = dense.qubit_count; + match gate { + &Gate::Single { op, qubit } => { + dense.apply1(qubit, gate_matrix(op)); + phased.left_mul(op, &[qubit]); + } + &Gate::Two { op, first, second } => { + match op { + ControlledX => dense.apply_cx(first, second), + ControlledZ => dense.apply_cz(first, second), + Swap => dense.apply_swap(first, second), + PrepareBell => { + dense.apply1(first, gate_matrix(Hadamard)); + dense.apply_cx(first, second); + } + _ => unreachable!("Gate::Two only carries two-qubit ops"), + } + phased.left_mul(op, &[first, second]); + } + Gate::Pauli(pauli) => { + let (x_bits, z_bits, phase) = pauli_arrays(pauli, qubit_count); + dense.apply_pauli(&x_bits, &z_bits, phase); + phased.left_mul_pauli(pauli); + } + Gate::PauliExp(pauli) => { + let (x_bits, z_bits, phase) = pauli_arrays(pauli, qubit_count); + dense.apply_pauli_exp(&x_bits, &z_bits, phase); + phased.left_mul_pauli_exp(pauli); + } + } +} + +proptest! { + #![proptest_config(ProptestConfig::with_cases(400))] + #[test] + fn phased_clifford_tracks_dense_statevector(gates in vec(gate_strategy(4), 0..40)) { + let qubit_count = 4; + let mut dense = Dense::zero(qubit_count); + let mut phased = PhasedCliffordUnitary::identity(qubit_count); + for gate in &gates { + apply(gate, &mut dense, &mut phased); + } + prop_assert!(close(&statevector(&phased), &dense.amp), "diverged on {gates:?}"); + } +} diff --git a/pauliverse/Cargo.toml b/pauliverse/Cargo.toml index c2a73cae..34463582 100644 --- a/pauliverse/Cargo.toml +++ b/pauliverse/Cargo.toml @@ -19,6 +19,7 @@ derive_more = { version = "2.0.1", features = [ [dev-dependencies] criterion = { version = "0.5", features = ["html_reports"] } proptest = "1.0" +dense-oracle = { path = "../test-utils/dense-oracle" } [[bench]] name = "simulation_benchmark" diff --git a/pauliverse/README.md b/pauliverse/README.md index 001b6ede..5c728b65 100644 --- a/pauliverse/README.md +++ b/pauliverse/README.md @@ -10,10 +10,11 @@ pauliverse provides multiple stabilizer simulation implementations optimized for - **`OutcomeSpecificSimulation`**: Traditional stabilizer simulation that draws random measurement outcomes as needed. Best for Monte Carlo sampling when the number of shots is much smaller than the number of random measurements. - **`OutcomeCompleteSimulation`**: Tracks all 2^n_random outcome branches simultaneously. Best for analyzing entire circuits, or when shots >> 2^n_random. +- **`PhasedOutcomeCompleteSimulation`**: Like `OutcomeCompleteSimulation`, but also tracks the exact global phase of the encoded state. Best for exact equality checking of non-stabilizer circuits. - **`OutcomeFreeSimulation`**: Tracks stabilizer modulo measurement outcomes. Best for circuit verification and logical operator analysis. - **`FaultySimulation`**: Extends OutcomeCompleteSimulation with frame-based noise propagation. Best for estimating logical error rates under Pauli noise models. -All simulators support the full Clifford group and Pauli measurements. Based on algorithms from [arXiv:2309.08676](https://arxiv.org/abs/2309.08676). +All simulators support the full Clifford group and Pauli measurements. Based on algorithms from [arXiv:2309.08676](https://arxiv.org/abs/2309.08676), with exact global-phase tracking from [arXiv:2603.24717](https://arxiv.org/abs/2603.24717). ## Installation @@ -91,6 +92,7 @@ outcome = sim.measure(pauli_z) |-----------|----------|---------------| | `OutcomeSpecificSimulation` | Monte Carlo sampling with many shots | Minimal overhead per shot, simple API | | `OutcomeCompleteSimulation` | Exact distributions, circuit verification | Simulates once, sample many times efficiently | +| `PhasedOutcomeCompleteSimulation` | Exact equality checking of non-stabilizer circuits | Tracks the exact global phase, not just up to phase | | `OutcomeFreeSimulation` | State verification, logical operators | Tracks stabilizers without measurement records | | `FaultySimulation` | Logical error rates, decoder testing | Frame-based noise propagation | @@ -134,6 +136,7 @@ Key resources: - [Simulation trait](src/lib.rs) - 40+ methods for gates, measurements, and state queries - [OutcomeSpecificSimulation](src/outcome_specific_simulation.rs) - Traditional simulation with random outcomes - [OutcomeCompleteSimulation](src/outcome_complete_simulation.rs) - All-branches simulation for exact analysis +- [PhasedOutcomeCompleteSimulation](src/phased_outcome_complete_simulation.rs) - All-branches simulation tracking the exact global phase - [FaultySimulation](src/faulty_simulation.rs) - Noisy simulation with frame propagation ## Contributing diff --git a/pauliverse/src/action.rs b/pauliverse/src/action.rs index e5acb6c6..63bfb517 100644 --- a/pauliverse/src/action.rs +++ b/pauliverse/src/action.rs @@ -1,7 +1,7 @@ use std::fmt::Debug; use crate::{ - OutcomeCompleteSimulation, Simulation, + OutcomeCompleteSimulation, PhasedOutcomeCompleteSimulation, Simulation, circuit::{Circuit, SimulationError}, }; use binar::{AffineMap, BitMatrix, BitVec, Bitwise, BitwiseMut, IndexSet}; @@ -84,6 +84,17 @@ pub enum ActionsInequivalenceReason { ChoiState, /// See [`CircuitAction::signed_choi_state_stabilizers`] for details. ChoiStateSigns, + /// The relative `ζ₈` phases between branches of the Choi state differ. + /// Only produced by [`PhasedCircuitAction`]; see its documentation for details. + RelativePhase, + /// The two phased actions have different numbers of symbolic-angle (virtual) random bits, so no + /// one-to-one correspondence between their symbolic rotations exists. + /// Only produced by [`PhasedCircuitAction`]; see its documentation for details. + SymbolicAngleCount, + /// A supplied outcome remapping would affinely mix symbolic-angle (virtual) random bits, either + /// with one another or with true (measurement) random bits, which does not correspond to any + /// operator equality. Only produced by [`PhasedCircuitAction::is_equivalent_with_map`]. + SymbolicAngleMixed, } /// [`Circuit`]s in pauliverse include fixed number of qubits and do not have prepare and destroy instructions. @@ -99,22 +110,96 @@ pub fn action_of( input_qubits: &[QubitId], output_qubits: &[QubitId], ) -> Result { + build_action::(circuit, input_qubits, output_qubits).map(|(action, _)| action) +} + +/// Stabilizer simulators that expose the encoder data required to compute a [`CircuitAction`]. +/// +/// The method names differ from the inherent accessors of the same purpose to avoid shadowing them +/// inside the forwarding implementations. +trait ActionSimulation: Simulation { + fn encoder(&self) -> CliffordUnitary; + fn signs(&self) -> BitMatrix; + fn random_indicator(&self) -> &[bool]; + fn outcomes(&self) -> BitMatrix; + fn outcome_offset(&self) -> BitVec; +} + +impl ActionSimulation for OutcomeCompleteSimulation { + fn encoder(&self) -> CliffordUnitary { + self.state_encoder() + } + fn signs(&self) -> BitMatrix { + self.sign_matrix() + } + fn random_indicator(&self) -> &[bool] { + self.random_outcome_indicator() + } + fn outcomes(&self) -> BitMatrix { + self.outcome_matrix() + } + fn outcome_offset(&self) -> BitVec { + self.outcome_shift() + } +} + +impl ActionSimulation for PhasedOutcomeCompleteSimulation { + fn encoder(&self) -> CliffordUnitary { + self.state_encoder() + } + fn signs(&self) -> BitMatrix { + self.sign_matrix() + } + fn random_indicator(&self) -> &[bool] { + self.random_outcome_indicator() + } + fn outcomes(&self) -> BitMatrix { + self.outcome_matrix() + } + fn outcome_offset(&self) -> BitVec { + self.outcome_shift() + } +} + +/// Computes a [`CircuitAction`] using simulator `S`, returning both the action and the consumed +/// simulator so that phase-aware callers can additionally read out its phase data. +fn build_action( + circuit: &Circuit, + input_qubits: &[QubitId], + output_qubits: &[QubitId], +) -> Result<(CircuitAction, S), ActionError> { let qubit_count = circuit .qubit_count() .max(input_qubits.iter().max().map_or(0, |&q| q + 1)) .max(output_qubits.iter().max().map_or(0, |&q| q + 1)); let reference_qubits: Vec = (qubit_count..qubit_count + input_qubits.len()).collect(); let outcome_count = circuit.outcome_count(); - let mut simulation = - OutcomeCompleteSimulation::with_capacity(qubit_count + input_qubits.len(), outcome_count, outcome_count); + let mut simulation = S::with_capacity(qubit_count + input_qubits.len(), outcome_count, outcome_count); for (input_qubit, reference_qubit) in input_qubits.iter().zip(reference_qubits.iter()) { simulation.unitary_op(paulimer::UnitaryOp::PrepareBell, &[*input_qubit, *reference_qubit]); } circuit.simulate(&mut simulation)?; - let sign_matrix = simulation.sign_matrix(); - let state_encoder = simulation.state_encoder(); + let action = action_from_simulation(&simulation, input_qubits, output_qubits, &reference_qubits, qubit_count)?; + Ok((action, simulation)) +} + +/// Canonicalizes the Choi state recorded in `simulation` into a [`CircuitAction`]. +/// +/// This is the post-simulation core shared by [`build_action`] (which prepares the Bell pairs and +/// replays a [`Circuit`]) and [`phased_action_from_simulation`] (which canonicalizes a Choi state the +/// caller has already prepared). The caller is responsible for having entangled `input_qubits[k]` +/// with `reference_qubits[k]` via a Bell pair before applying the circuit. +fn action_from_simulation( + simulation: &S, + input_qubits: &[QubitId], + output_qubits: &[QubitId], + reference_qubits: &[QubitId], + qubit_count: usize, +) -> Result { + let sign_matrix = simulation.signs(); + let state_encoder = simulation.encoder(); let auxiliary_qubits: Vec = output_qubits .iter() @@ -132,7 +217,7 @@ pub fn action_of( }); } - let observables = GeneratorsWithSigns::from_restriction(&state_encoder, &sign_matrix, &reference_qubits, true); + let observables = GeneratorsWithSigns::from_restriction(&state_encoder, &sign_matrix, reference_qubits, true); let stabilizers = GeneratorsWithSigns::from_restriction(&state_encoder, &sign_matrix, output_qubits, false); let choi_state_stabilizers = GeneratorsWithSigns::from_restriction( &state_encoder, @@ -145,13 +230,13 @@ pub fn action_of( false, ); - let indicators = simulation.random_outcome_indicator(); + let indicators = simulation.random_indicator(); let random_bit_map_matrix = random_bit_map_matrix(indicators); - let random_bit_map_shift = &random_bit_map_matrix * &simulation.outcome_shift().as_view(); + let random_bit_map_shift = &random_bit_map_matrix * &simulation.outcome_offset().as_view(); let outcome_to_random_bit_map = AffineMap::affine(random_bit_map_matrix.clone(), random_bit_map_shift.clone()); - let outcomes_from_random = AffineMap::affine(simulation.outcome_matrix(), simulation.outcome_shift().clone()); + let outcomes_from_random = AffineMap::affine(simulation.outcomes(), simulation.outcome_offset().clone()); - let action = CircuitAction { + Ok(CircuitAction { observables, stabilizers, choi_state_stabilizers, @@ -159,8 +244,7 @@ pub fn action_of( random_from_outcomes: outcome_to_random_bit_map, outcomes_from_random, input_qubit_ids: input_qubits.to_vec(), - }; - Ok(action) + }) } impl CircuitAction { @@ -325,10 +409,435 @@ impl CircuitAction { } } +/// The exact-global-phase analog of [`CircuitAction`], computed with a +/// [`PhasedOutcomeCompleteSimulation`] so that the **relative `ζ₈` phases between branches** of the +/// circuit's Choi state are retained in addition to the phaseless stabilizer data. +/// +/// A [`CircuitAction`] determines the Choi state only up to phase, so it cannot distinguish circuits +/// that act identically on the Pauli group but differ by branch-dependent phases — for example +/// `e^{iα Z}` and `e^{-iα Z}`, whose conditioned Paulis `+Z` and `-Z` share a symplectic action. +/// [`PhasedCircuitAction`] additionally compares the per-branch phase function +/// `φ(r) = i^⟨p, r⟩ (-1)^⟨B r + s, r⟩`, capturing exactly that information. +/// +/// The comparison is *up to a single global phase* common to all branches: the encoder's absolute +/// phase is not exposed, so two Choi states that differ only by an overall scalar are reported as +/// equivalent. Pinning down that absolute phase as well requires the auxiliary-qubit separation of +/// §4.5 of [arXiv:2603.24717](https://arxiv.org/abs/2603.24717), a planned follow-up. +#[derive(Debug, Clone, PartialEq)] +pub struct PhasedCircuitAction { + action: CircuitAction, + phase: PhaseData, + /// Indicator over the inner random bits: `true` where the bit is a symbolic rotation angle (a + /// "virtual" random bit allocated via [`Simulation::allocate_symbolic_angle`]) rather than a + /// genuine measurement-derived random bit. + symbolic_angles: BitVec, +} + +/// Computes a [`PhasedCircuitAction`] for `circuit` with the given input and output qubits. +/// +/// Behaves exactly like [`action_of`] but uses a [`PhasedOutcomeCompleteSimulation`], additionally +/// recording the branch phase function of the circuit's Choi state. +/// +/// # Errors +/// +/// Returns [`ActionError`] if action calculation fails. +pub fn phased_action_of( + circuit: &Circuit, + input_qubits: &[QubitId], + output_qubits: &[QubitId], +) -> Result { + let (action, simulation) = build_action::(circuit, input_qubits, output_qubits)?; + Ok(phased_action(action, &simulation)) +} + +/// Computes a [`PhasedCircuitAction`] directly from a [`PhasedOutcomeCompleteSimulation`] whose Choi +/// state the caller has already prepared. +/// +/// This is the simulator-native counterpart of [`phased_action_of`], matching the convention used by +/// the Python bindings where the simulator itself records the circuit. The caller must, before +/// applying the circuit, have entangled each `input_qubits[k]` with a reference qubit via +/// `UnitaryOp::PrepareBell`, following the same layout as [`phased_action_of`]: the reference qubit +/// for `input_qubits[k]` is `system_qubit_count + k`, where `system_qubit_count` is one past the +/// largest index appearing in `input_qubits` or `output_qubits`. +/// +/// # Errors +/// +/// Returns [`ActionError::AuxiliaryQubitsEntangled`] if the non-output system qubits remain +/// entangled with the rest of the state. +pub fn phased_action_from_simulation( + simulation: &PhasedOutcomeCompleteSimulation, + input_qubits: &[QubitId], + output_qubits: &[QubitId], +) -> Result { + let system_qubit_count = input_qubits + .iter() + .chain(output_qubits.iter()) + .copied() + .max() + .map_or(0, |qubit| qubit + 1); + let reference_qubits: Vec = (system_qubit_count..system_qubit_count + input_qubits.len()).collect(); + let action = action_from_simulation( + simulation, + input_qubits, + output_qubits, + &reference_qubits, + system_qubit_count, + )?; + Ok(phased_action(action, simulation)) +} + +/// Assembles a [`PhasedCircuitAction`] from a computed `action` and the `simulation` that recorded +/// the branch phase function. +fn phased_action(action: CircuitAction, simulation: &PhasedOutcomeCompleteSimulation) -> PhasedCircuitAction { + let symbolic_angles: BitVec = simulation.symbolic_angle_indicator().iter().copied().collect(); + PhasedCircuitAction { + action, + phase: PhaseData::from_simulation(simulation), + symbolic_angles, + } +} + +impl PhasedCircuitAction { + /// The underlying phaseless [`CircuitAction`]. + #[must_use] + pub fn action(&self) -> &CircuitAction { + &self.action + } + + /// Canonical choi state stabilizers; see [`CircuitAction::choi_state_stabilizers`]. + pub fn choi_state_stabilizers(&self) -> &[SparsePauli] { + self.action.choi_state_stabilizers() + } + + /// Returns `Ok(())` if the phaseless actions are equivalent up to signs, otherwise the reasons. + /// + /// This ignores phase entirely; use [`Self::is_equivalent_with_map`] to additionally compare the + /// relative branch phases. + /// + /// # Errors + /// + /// Returns a list of [`ActionsInequivalenceReason`] if the phaseless actions differ. + pub fn is_equivalent_up_to_signs( + &self, + other: &PhasedCircuitAction, + ) -> Result<(), Vec> { + self.action.is_equivalent_up_to_signs(&other.action) + } + + /// Verifies that two phased actions implement the same operator on every input, enforcing the + /// **virtual/true random-bit distinction**: symbolic-angle (virtual) random bits must correspond + /// *one to one* between the two actions, while true (measurement-derived) random bits may be + /// marginalized. + /// + /// A symbolic rotation `e^{iα P}` is modelled by conditioning `P` on a bit allocated via + /// [`Simulation::allocate_symbolic_angle`]. Two encodings of the same parameterised circuit are + /// equivalent only when their angle bits match up identically — angle `α_k` of one maps to angle + /// `α_k` of the other, in allocation order, with no affine mixing. True random bits (allocated + /// via [`Simulation::allocate_random_bit`] or produced by a genuine measurement) carry no such + /// constraint: surplus true bits present in only one action are projected out, matching the way + /// the phaseless [`CircuitAction::is_equivalent_with_map`] marginalizes measurement randomness. + /// This is exactly what makes measurement-based "ejection" gadgets compare equal to the operation + /// they implement directly. + /// + /// The two actions must have the same number of symbolic angles (otherwise + /// [`ActionsInequivalenceReason::SymbolicAngleCount`] is returned). When both actions also share + /// the same true random bits and those bits must be related non-trivially, use + /// [`Self::is_equivalent_with_map`] with an explicit correspondence. + /// + /// # Errors + /// + /// Returns a list of [`ActionsInequivalenceReason`] if the actions differ. + pub fn is_equivalent(&self, other: &PhasedCircuitAction) -> Result<(), Vec> { + let map = self.provenance_random_map(other).map_err(|reason| vec![reason])?; + self.check_with_random_map(other, &map) + } + + /// Check if two phased actions are equivalent (up to a single global phase) when outcomes are + /// remapped, comparing both the [`CircuitAction`] data and the relative branch phases. + /// + /// The outcome remapping `self_outcomes_from_other_outcomes` follows the same convention as + /// [`CircuitAction::is_equivalent_with_map`]: outcomes of `self` equal `A(o_other)`. When the map + /// is `None`, the zero map is used, as is common for circuits with unitary action. + /// + /// This is the lower-level escape hatch behind [`Self::is_equivalent`]. The supplied map may + /// affinely remap *true* random bits, but **symbolic-angle (virtual) random bits must be mapped + /// one to one** in allocation order. Any map whose induced random-bit correspondence affinely + /// combines angle bits, or mixes them with true random bits, is rejected with + /// [`ActionsInequivalenceReason::SymbolicAngleMixed`]; prefer [`Self::is_equivalent`] unless you + /// specifically need to relabel true random bits. + /// + /// # Errors + /// + /// Returns a list of [`ActionsInequivalenceReason`] if the actions differ; the additional + /// [`ActionsInequivalenceReason::RelativePhase`] is returned when only the branch phases differ. + pub fn is_equivalent_with_map( + &self, + other: &PhasedCircuitAction, + self_outcomes_from_other_outcomes: Option<&AffineMap>, + ) -> Result<(), Vec> { + let zero = zero_map(&self.action, &other.action); + let outcome_map = self_outcomes_from_other_outcomes.unwrap_or(&zero); + let self_outcomes_from_other_random = outcome_map.dot(&other.action.outcomes_from_random); + let self_random_from_other_random = self.action.random_from_outcomes.dot(&self_outcomes_from_other_random); + + if !self.angle_correspondence_is_clean(other, &self_random_from_other_random) { + return Err(vec![ActionsInequivalenceReason::SymbolicAngleMixed]); + } + self.check_with_random_map(other, &self_random_from_other_random) + } + + /// Builds the random-bit correspondence used by [`Self::is_equivalent`]: identity (in allocation + /// order) on the symbolic-angle bits, identity on the true bits shared by both actions, and a + /// projection to zero of any surplus true bits present only in `other`. + fn provenance_random_map(&self, other: &PhasedCircuitAction) -> Result { + let self_angles: Vec = self.symbolic_angles.support().collect(); + let other_angles: Vec = other.symbolic_angles.support().collect(); + if self_angles.len() != other_angles.len() { + return Err(ActionsInequivalenceReason::SymbolicAngleCount); + } + let self_random = self.symbolic_angles.len(); + let other_random = other.symbolic_angles.len(); + let self_trues = (0..self_random).filter(|&index| !self.symbolic_angles.index(index)); + let other_trues: Vec = (0..other_random) + .filter(|&index| !other.symbolic_angles.index(index)) + .collect(); + + let mut matrix = BitMatrix::zeros(self_random, other_random); + for (&self_bit, &other_bit) in self_angles.iter().zip(other_angles.iter()) { + matrix.set((self_bit, other_bit), true); + } + for (self_bit, &other_bit) in self_trues.zip(other_trues.iter()) { + matrix.set((self_bit, other_bit), true); + } + Ok(AffineMap::linear(matrix)) + } + + /// Runs the count, sign and relative-phase comparisons under a given random-bit correspondence + /// `self_random_from_other_random` (branch `r` of `other` corresponds to branch + /// `self_random_from_other_random(r)` of `self`). + fn check_with_random_map( + &self, + other: &PhasedCircuitAction, + self_random_from_other_random: &AffineMap, + ) -> Result<(), Vec> { + self.action.is_equivalent_up_to_signs(&other.action)?; + + let mut reasons = Vec::new(); + if self + .action + .observables + .is_equivalent_with_map(&other.action.observables, self_random_from_other_random) + { + reasons.push(ActionsInequivalenceReason::ObservablesSigns); + } + if self + .action + .stabilizers + .is_equivalent_with_map(&other.action.stabilizers, self_random_from_other_random) + { + reasons.push(ActionsInequivalenceReason::StabilizersSigns); + } + if self + .action + .choi_state_stabilizers + .is_equivalent_with_map(&other.action.choi_state_stabilizers, self_random_from_other_random) + { + reasons.push(ActionsInequivalenceReason::ChoiStateSigns); + } + if !self.relative_phase_matches(other) { + reasons.push(ActionsInequivalenceReason::RelativePhase); + } + if reasons.is_empty() { Ok(()) } else { Err(reasons) } + } + + /// Guards against an outcome remapping that does not respect the virtual/true distinction. + /// + /// Returns `true` iff the induced random-bit correspondence maps the symbolic-angle bits of + /// `other` one to one onto those of `self` (in allocation order) with no leakage: each angle bit + /// of `other` maps exactly to the matching angle bit of `self`, and no angle bit of `self` is + /// driven by a true (non-angle) bit of `other`. + fn angle_correspondence_is_clean( + &self, + other: &PhasedCircuitAction, + self_random_from_other_random: &AffineMap, + ) -> bool { + let self_angles: Vec = self.symbolic_angles.support().collect(); + let other_angles: Vec = other.symbolic_angles.support().collect(); + if self_angles.len() != other_angles.len() { + return false; + } + let matrix = self_random_from_other_random.matrix(); + if self_angles + .iter() + .any(|&self_bit| self_random_from_other_random.shift().index(self_bit)) + { + return false; + } + for (&self_angle, &other_angle) in self_angles.iter().zip(other_angles.iter()) { + for self_bit in 0..matrix.row_count() { + let expected = self_bit == self_angle; + if matrix.get((self_bit, other_angle)) != expected { + return false; + } + } + } + for &self_angle in &self_angles { + for other_bit in 0..matrix.column_count() { + if !other_angles.contains(&other_bit) && matrix.get((self_angle, other_bit)) { + return false; + } + } + } + true + } + + /// Checks that the branch phase functions of `self` and `other` agree up to a global phase + /// **on the symbolic-angle (virtual) random bits only**. + /// + /// Symbolic-angle bits carry the coherent, observable relative phases of the modelled rotations + /// `e^{iα P}`, so they must match. True (measurement-derived) random bits label incoherent, + /// traced-out measurement branches whose per-branch global phase is physically unobservable, so + /// the comparison ignores them entirely (it sets every true bit to zero). This is what makes the + /// phased equivalence reduce *exactly* to the phaseless [`CircuitAction`] equivalence when no + /// symbolic angles are present, and what lets measurement-based "ejection" gadgets — whose + /// corrected ancilla branches differ only by an unobservable per-branch phase — compare equal to + /// the operation they implement directly. + /// + /// The angle bits of `self` and `other` correspond one to one in allocation order. The phase + /// function `φ(r) = 2⟨p, r⟩ + 4⟨B r + s, r⟩ (mod 8)` restricted to the angle subspace is a + /// degree-≤2 polynomial, so it is fully determined by its values on the zero vector, the angle + /// unit vectors, and their pairwise sums. Equality up to a global phase is therefore equivalent + /// to equality of the linear coefficients `φ(e_i) − φ(0)` and the quadratic coefficients + /// `φ(e_i + e_j) − φ(e_i) − φ(e_j) + φ(0)` (ignoring the constant `φ(0)`, i.e. the global phase). + fn relative_phase_matches(&self, other: &PhasedCircuitAction) -> bool { + let self_angles: Vec = self.symbolic_angles.support().collect(); + let other_angles: Vec = other.symbolic_angles.support().collect(); + if self_angles.len() != other_angles.len() { + return false; + } + let angle_count = self_angles.len(); + let self_dimension = self.phase.random_count(); + let other_dimension = other.phase.random_count(); + + let phase_self = |angles: &[usize]| { + let indices: Vec = angles.iter().map(|&order| self_angles[order]).collect(); + self.phase.phase_exponent(&unit_vector(self_dimension, &indices)) + }; + let phase_other = |angles: &[usize]| { + let indices: Vec = angles.iter().map(|&order| other_angles[order]).collect(); + other.phase.phase_exponent(&unit_vector(other_dimension, &indices)) + }; + + let constant_self = i32::from(phase_self(&[])); + let constant_other = i32::from(phase_other(&[])); + + let mut linear_self = vec![0i32; angle_count]; + let mut linear_other = vec![0i32; angle_count]; + for order in 0..angle_count { + linear_self[order] = (i32::from(phase_self(&[order])) - constant_self).rem_euclid(8); + linear_other[order] = (i32::from(phase_other(&[order])) - constant_other).rem_euclid(8); + } + if linear_self != linear_other { + return false; + } + + for first in 0..angle_count { + for second in (first + 1)..angle_count { + let quadratic_self = (i32::from(phase_self(&[first, second])) + - constant_self + - linear_self[first] + - linear_self[second]) + .rem_euclid(8); + let quadratic_other = (i32::from(phase_other(&[first, second])) + - constant_other + - linear_other[first] + - linear_other[second]) + .rem_euclid(8); + if quadratic_self != quadratic_other { + return false; + } + } + } + true + } +} + // ================================================================================================ // Private Types // ================================================================================================ +/// Branch phase function of a Choi state, indexed by the inner random bits. +/// +/// The `ζ₈` phase of branch `r` is `ζ₈^φ(r)` with `φ(r) = 2⟨p, r⟩ + 4⟨B r + s, r⟩ (mod 8)`, matching +/// [`PhasedOutcomeCompleteSimulation::output_phase_exponent`]. +#[derive(Debug, Clone, PartialEq)] +pub(crate) struct PhaseData { + /// `p`: linear `i` phase. + linear_i: BitVec, + /// `s`: linear `-1` phase. + linear_sign: BitVec, + /// `B`: quadratic `-1` phase. + quadratic: BitMatrix, +} + +impl PhaseData { + /// Extracts the branch phase function recorded by `simulation`. + pub(crate) fn from_simulation(simulation: &PhasedOutcomeCompleteSimulation) -> Self { + PhaseData { + linear_i: simulation.linear_i_phase(), + linear_sign: simulation.linear_sign_phase(), + quadratic: simulation.quadratic_phase_matrix(), + } + } + + fn random_count(&self) -> usize { + self.linear_i.len() + } + + /// The `ζ₈` exponent `φ(r) = 2⟨p, r⟩ + 4⟨B r + s, r⟩ (mod 8)` for the branch `random_bits`. + pub(crate) fn phase_exponent(&self, random_bits: &BitVec) -> u8 { + phase_form_exponent( + self.random_count(), + |index| random_bits.index(index), + |index| self.linear_i.index(index), + |index| self.linear_sign.index(index), + |row, column| self.quadratic.get((row, column)), + ) + } +} + +/// Evaluates the `ζ₈ = e^{iπ/4}` exponent of the F₂ phase form `i^⟨p, r⟩ (-1)^⟨B r + s, r⟩` for a +/// random-bit assignment `r`. +/// +/// The coefficients are read through accessor closures so the phased simulator and its lowered +/// [`crate::action`] `PhaseData` — which store `p`, `s`, `B` and `r` in different (aligned vs. +/// unaligned) representations — share a single implementation. `random_bit`, `linear_i` (`p`) and +/// `linear_sign` (`s`) are indexed by column and `quadratic` reads `B[(row, column)]`, all over +/// `0..random_count`. +pub(crate) fn phase_form_exponent( + random_count: usize, + random_bit: impl Fn(usize) -> bool, + linear_i: impl Fn(usize) -> bool, + linear_sign: impl Fn(usize) -> bool, + quadratic: impl Fn(usize, usize) -> bool, +) -> u8 { + let mut linear_i_parity = false; + let mut sign = false; + for column in 0..random_count { + if !random_bit(column) { + continue; + } + linear_i_parity ^= linear_i(column); + sign ^= linear_sign(column); + for row in 0..random_count { + if random_bit(row) && quadratic(row, column) { + sign = !sign; + } + } + } + (2 * u8::from(linear_i_parity) + 4 * u8::from(sign)) % 8 +} + #[derive(Debug, Clone, PartialEq)] struct GeneratorsWithSigns { /// Canonical choice of generators, with canonical signs @@ -421,3 +930,12 @@ fn adjust_phase_to_canonical(pauli: &mut SparsePauli) -> bool { fn zero_map(to: &CircuitAction, from: &CircuitAction) -> AffineMap { AffineMap::zero(from.outcome_count(), to.outcome_count()) } + +/// Returns the length-`dimension` bit vector with the bits in `set_indices` set to one. +fn unit_vector(dimension: usize, set_indices: &[usize]) -> BitVec { + let mut vector = BitVec::zeros(dimension); + for &index in set_indices { + vector.assign_index(index, true); + } + vector +} diff --git a/pauliverse/src/circuit.rs b/pauliverse/src/circuit.rs index 31f46e17..29c82dc1 100644 --- a/pauliverse/src/circuit.rs +++ b/pauliverse/src/circuit.rs @@ -39,6 +39,9 @@ pub(crate) enum Instruction { }, AllocateRandomBit { outcome_id: OutcomeId, + /// Whether this bit is a symbolic rotation angle (a "virtual" random bit) rather than a + /// genuine random bit. See [`Simulation::allocate_symbolic_angle`]. + symbolic_angle: bool, }, ConditionalPauli { pauli: SparsePauli, @@ -203,8 +206,15 @@ impl Circuit { }); } } - Instruction::AllocateRandomBit { outcome_id } => { - let sim_outcome_id = simulator.allocate_random_bit(); + Instruction::AllocateRandomBit { + outcome_id, + symbolic_angle, + } => { + let sim_outcome_id = if *symbolic_angle { + simulator.allocate_symbolic_angle() + } else { + simulator.allocate_random_bit() + }; if *outcome_id != sim_outcome_id { return Err(SimulationError::InvalidInstructionOutcomeId { expected: *outcome_id, @@ -348,7 +358,7 @@ impl CircuitBuilder { /// Push a raw instruction to the circuit. pub(crate) fn push(&mut self, instruction: Instruction) { match instruction { - Instruction::Measure { outcome_id, .. } | Instruction::AllocateRandomBit { outcome_id } => { + Instruction::Measure { outcome_id, .. } | Instruction::AllocateRandomBit { outcome_id, .. } => { assert_eq!( outcome_id, self.outcome_count, "Instruction outcome_id {outcome_id} does not match expected outcome_count" @@ -365,7 +375,20 @@ impl Simulation for CircuitBuilder { fn allocate_random_bit(&mut self) -> OutcomeId { let outcome_id = self.outcome_count; self.outcome_count += 1; - self.circuit.push(Instruction::AllocateRandomBit { outcome_id }); + self.circuit.push(Instruction::AllocateRandomBit { + outcome_id, + symbolic_angle: false, + }); + outcome_id + } + + fn allocate_symbolic_angle(&mut self) -> OutcomeId { + let outcome_id = self.outcome_count; + self.outcome_count += 1; + self.circuit.push(Instruction::AllocateRandomBit { + outcome_id, + symbolic_angle: true, + }); outcome_id } @@ -604,7 +627,10 @@ mod tests { _ => { let outcome_id = *outcome_counter; *outcome_counter += 1; - Instruction::AllocateRandomBit { outcome_id } + Instruction::AllocateRandomBit { + outcome_id, + symbolic_angle: false, + } } } } @@ -792,7 +818,10 @@ mod tests { #[test] fn allocate_random_bit_has_no_faults() { let mut circuit = Circuit::new(); - circuit.push(Instruction::AllocateRandomBit { outcome_id: 0 }); + circuit.push(Instruction::AllocateRandomBit { + outcome_id: 0, + symbolic_angle: false, + }); assert_eq!(circuit.fault_count(), 0); assert_eq!(circuit.outcome_count(), 1); } diff --git a/pauliverse/src/lib.rs b/pauliverse/src/lib.rs index ec4f2274..3692340a 100644 --- a/pauliverse/src/lib.rs +++ b/pauliverse/src/lib.rs @@ -4,11 +4,12 @@ //! for different use cases in quantum computing and quantum error correction. //! //! These simulation algorithms are based on the framework described in -//! [arXiv:2309.08676](https://arxiv.org/abs/2309.08676). +//! [arXiv:2309.08676](https://arxiv.org/abs/2309.08676), extended with exact +//! global-phase tracking from [arXiv:2603.24717](https://arxiv.org/abs/2603.24717). //! //! # Overview //! -//! This crate offers four simulation modes: +//! This crate offers five simulation modes: //! //! - **[`OutcomeSpecificSimulation`]**: Traditional simulation with random (or caller-supplied) measurement outcomes. //! Best for Monte Carlo sampling and estimating error rates. @@ -16,6 +17,11 @@ //! - **[`OutcomeCompleteSimulation`]**: Tracks all possible measurement outcomes simultaneously. //! Achieves asymptotic speedup when enumerating outcomes. //! +//! - **[`PhasedOutcomeCompleteSimulation`]**: Like [`OutcomeCompleteSimulation`], but additionally +//! tracks the *exact* global phase of the encoded state (Algorithm 4.2 of +//! [arXiv:2603.24717](https://arxiv.org/abs/2603.24717)). Enables exact equality checking of +//! non-stabilizer circuits. +//! //! - **[`OutcomeFreeSimulation`]**: Simulation without tracking specific outcomes. //! Minimal overhead when you only care about stabilizer state evolution up to global phase. //! @@ -51,6 +57,7 @@ //! |-----------|----------|---------------| //! | [`OutcomeSpecificSimulation`] | Monte Carlo sampling with few shots | Single concrete execution path | //! | [`OutcomeCompleteSimulation`] | Whole-circuit analysis, enumerating outcomes | Avoids re-simulating for each outcome sample | +//! | [`PhasedOutcomeCompleteSimulation`] | Exact equality checking of non-stabilizer circuits | Tracks the exact global phase, not just up to phase | //! | [`OutcomeFreeSimulation`] | Stabilizer queries without outcomes | Minimal overhead, no outcome tracking | //! | [`FaultySimulation`] | Noisy simulation with error correction | Efficient frame-based noise propagation | //! @@ -97,6 +104,7 @@ pub mod outcome_complete_simulation; pub mod outcome_free_simulation; pub mod outcome_specific_simulation; +pub mod phased_outcome_complete_simulation; pub mod sampling; #[cfg(test)] pub(crate) mod statistical_testing; @@ -108,6 +116,7 @@ pub use noise::{OutcomeCondition, PauliDistribution, PauliFault}; pub use outcome_complete_simulation::OutcomeCompleteSimulation; pub use outcome_free_simulation::OutcomeFreeSimulation; pub use outcome_specific_simulation::OutcomeSpecificSimulation; +pub use phased_outcome_complete_simulation::PhasedOutcomeCompleteSimulation; type Pauli = paulimer::pauli::SparsePauli; type Unitary = paulimer::clifford::CliffordUnitary; @@ -129,6 +138,39 @@ pub trait Simulation: Default { /// Returns the outcome ID for the newly allocated outcome. fn allocate_random_bit(&mut self) -> OutcomeId; + /// Allocate a new outcome representing a *symbolic rotation angle*. + /// + /// A symbolic angle is the parameter of an exponential `exp(iαP)`, realised by applying `P` + /// conditioned on the returned outcome (e.g. `conditional_pauli(P, &[angle], true)`). It behaves + /// like a random bit during simulation, but it carries different *provenance*: symbolic angles + /// model the (unknown) continuous parameters of a circuit, whereas [`Self::allocate_random_bit`] + /// models genuine measurement randomness. + /// + /// Simulators that track this distinction (such as the phased outcome-complete simulator) use it + /// to require that symbolic angles correspond one-to-one between circuits being compared, while + /// genuine random bits may be affinely remapped. The default implementation simply defers to + /// [`Self::allocate_random_bit`], so simulators that do not track provenance are unaffected. + /// + /// Returns the outcome ID for the newly allocated symbolic angle. + fn allocate_symbolic_angle(&mut self) -> OutcomeId { + self.allocate_random_bit() + } + + /// Apply a symbolic Pauli rotation `exp(iα P)`, parameterised by the symbolic `angle` `α`. + /// + /// This is the high-level way to add a parameterised rotation to a circuit: allocate the angle + /// with [`Self::allocate_symbolic_angle`], then call this with the Pauli `P`. Prefer it over + /// manually conditioning `P` on the angle bit via [`Self::conditional_pauli`] — it states the + /// intent (a symbolic rotation) directly and keeps the angle's special provenance explicit. + /// + /// Because the comparison of two circuits matches symbolic angles one-to-one, the same `angle` + /// can parameterise several rotations to model a *shared* parameter `α`, and labelling the + /// rotations of two circuits with angles allocated in the same order is what makes them + /// correspond when the circuits are compared for equivalence. + fn symbolic_pauli_exp(&mut self, observable: &Pauli, angle: OutcomeId) { + self.conditional_pauli(observable, &[angle], true); + } + // ========== Unitary Operations ========== /// Apply a Clifford unitary to specified qubits. diff --git a/pauliverse/src/phased_outcome_complete_simulation.rs b/pauliverse/src/phased_outcome_complete_simulation.rs new file mode 100644 index 00000000..9ed43662 --- /dev/null +++ b/pauliverse/src/phased_outcome_complete_simulation.rs @@ -0,0 +1,662 @@ +use crate::Simulation; +use crate::action::phase_form_exponent; +use crate::outcome_complete_simulation::row_sum; +use crate::outcome_free_simulation::{max_pair_support, max_support}; +use binar::{BitMatrix, BitVec}; +use binar::{Bitwise, BitwiseMut, BitwisePair, BitwisePairMut, IndexSet, matrix::AlignedBitMatrix, vec::AlignedBitVec}; +use paulimer::clifford::{Clifford, CliffordUnitary, PhasedCliffordUnitary}; +use paulimer::pauli::{Pauli, PauliBits, PauliUnitary, anti_commutes_with, generic::PhaseExponent}; +use paulimer::pauli::{PauliBinaryOps, PauliMutable}; +use paulimer::{CLIFFORD_BIT_ALIGNMENT, UnitaryOp}; +use rand::RngExt; + +type SparsePauli = paulimer::pauli::SparsePauli; + +/// Asymptotically efficient stabilizer simulation tracking all outcomes *and* the exact global phase. +/// +/// This is the global-phase-tracking generalization of [`crate::OutcomeCompleteSimulation`]. Where +/// the latter implements Algorithm 5.3 of [KBP](https://arxiv.org/abs/2309.08676) and represents the +/// simulated state only *up to a global phase*, this simulator implements Algorithm 4.2 of the +/// [phased simulation paper](https://arxiv.org/abs/2603.24717) and tracks the global phase exactly. +/// +/// Exact phase tracking is what enables verification of circuits that contain *non-stabilizer* +/// resources such as symbolic single-qubit rotations: an equality `C₁ e^{iαZ} C₂|0⟩ = D₁ e^{iαZ} D₂|0⟩` +/// for all `α` reduces to a pair of exact stabilizer-state equalities, and exactness — not equality +/// up to a global phase — is precisely what makes the reduction valid. +/// +/// # Representation +/// +/// For a circuit with `n` output qubits and `n_M` outcomes the simulator maintains a vector +/// `q ∈ {1, 1/2}^{n_M}` of conditional outcome probabilities together with a *phased* Clifford +/// encoder `R` and `𝔽₂` data `A`, `B`, `M`, `v₀`, `p`, `s`. For a random-bit vector +/// `r ∈ {0,1}^{n_r}` (where `n_r` is the number of random outcomes) the outcome vector is +/// `v = v₀ + M r` and the exact output state is +/// +/// ```text +/// i^⟨p, r⟩ · (-1)^⟨B r + s, r⟩ · R |A r⟩. +/// ``` +/// +/// Here `⟨p, r⟩` is linear in `r` (an `i`-phase), while `⟨B r + s, r⟩` is a quadratic form in `r` +/// (a `±1`-phase). The encoder `R` additionally fixes the overall `ζ₈ = e^{iπ/4}` phase of `R|Ar⟩`. +/// +/// | code field | paper symbol | meaning | +/// |------------------------|--------------|-----------------------------------------------| +/// | `phased_clifford` | `R` | phased state encoder (exact global phase) | +/// | `sign_matrix` | `A` | random bits → computational-basis register | +/// | `quadratic_phase_matrix` | `B` | quadratic `-1` phase | +/// | `outcome_matrix` | `M` | random bits → outcome vector | +/// | `outcome_shift` | `v₀` | deterministic outcome shift | +/// | `linear_i_phase` | `p` | linear `i` phase | +/// | `linear_sign_phase` | `s` | linear `-1` phase | +/// | `random_outcome_indicator` | `q` | which outcomes are random (probability 1/2) | +/// +/// # Phase-resolved Clifford application +/// +/// A bare [`CliffordUnitary`] is *phaseless*: its symplectic tableau fixes the operator on the Pauli +/// group but not the overall `ζ₈` factor. Exact phase tracking therefore requires Cliffords to be +/// applied through phase-resolved entry points — the named gates of [`Simulation::unitary_op`], +/// [`Simulation::pauli`] and [`Simulation::pauli_exp`], plus [`Simulation::permute`], +/// [`Simulation::conditional_pauli`] and the measurement methods. Applying an opaque phaseless +/// [`CliffordUnitary`] via [`Simulation::clifford`] would leave the global phase undetermined, so +/// that method is unsupported here (it panics); see its documentation for the rationale. +/// +/// # Examples +/// +/// ``` +/// use pauliverse::{PhasedOutcomeCompleteSimulation, Simulation}; +/// use paulimer::UnitaryOp; +/// +/// let mut sim = PhasedOutcomeCompleteSimulation::new(2); +/// sim.unitary_op(UnitaryOp::Hadamard, &[0]); +/// sim.unitary_op(UnitaryOp::ControlledX, &[0, 1]); +/// +/// // The encoder now prepares a Bell pair with an exactly-known global phase. +/// assert_eq!(sim.random_outcome_count(), 0); +/// ``` +#[must_use] +pub struct PhasedOutcomeCompleteSimulation { + phased_clifford: PhasedCliffordUnitary, // R (phased encoder) + sign_matrix: AlignedBitMatrix, // A + quadratic_phase_matrix: AlignedBitMatrix, // B + outcome_matrix: AlignedBitMatrix, // M + outcome_shift: AlignedBitVec, // v_0 + linear_i_phase: AlignedBitVec, // p + linear_sign_phase: AlignedBitVec, // s + random_outcome_indicator: Vec, // vec(q), [j] is true iff vec(q)_j = 1/2 + symbolic_angle_indicator: Vec, // [k] is true iff random bit k is a symbolic rotation angle + random_bit_count: usize, + qubit_count: usize, +} + +impl std::fmt::Debug for PhasedOutcomeCompleteSimulation { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("PhasedOutcomeCompleteSimulation") + .field("state_encoder", &self.phased_clifford.clone().into_clifford()) + .field("sign_matrix", &self.sign_matrix()) + .field("quadratic_phase_matrix", &self.quadratic_phase_matrix()) + .field("outcome_matrix", &self.outcome_matrix()) + .field("outcome_shift", &self.outcome_shift().iter().collect::>()) + .field("linear_i_phase", &self.linear_i_phase().iter().collect::>()) + .field( + "linear_sign_phase", + &self.linear_sign_phase().iter().collect::>(), + ) + .field("random_outcome_indicator", &self.random_outcome_indicator) + .field("symbolic_angle_indicator", &self.symbolic_angle_indicator) + .field("random_bit_count", &self.random_bit_count) + .field("qubit_count", &self.qubit_count) + .finish() + } +} + +impl Default for PhasedOutcomeCompleteSimulation { + fn default() -> Self { + PhasedOutcomeCompleteSimulation::with_capacity(0, 0, 0) + } +} + +impl PhasedOutcomeCompleteSimulation { + /// Get the phaseless Clifford unitary encoding the current stabilizer state. + /// + /// This is the unitary `R` such that `R|0⟩` represents the stabilizer state, with the global + /// phase discarded. Use [`Self::phased_state_encoder`] to retain the exact global phase. + pub fn state_encoder(&self) -> CliffordUnitary { + self.phased_state_encoder().into_clifford() + } + + /// Get the phased Clifford unitary encoding the current stabilizer state with exact global phase. + pub fn phased_state_encoder(&self) -> PhasedCliffordUnitary { + let mut res = self.phased_clifford.clone(); + res.resize(self.qubit_count); + res + } + + /// Get the sign matrix `A` tracking how computational-basis registers depend on random bits. + /// + /// Returns a cache-aligned reference for efficiency. + pub fn aligned_sign_matrix(&self) -> &AlignedBitMatrix { + &self.sign_matrix + } + + /// Get a copy of the sign matrix `A` without alignment constraints. + pub fn sign_matrix(&self) -> BitMatrix { + BitMatrix::from_aligned(AlignedBitMatrix::from_row_iter( + self.sign_matrix.row_iterator(0..self.qubit_count()), + self.random_outcome_count(), + )) + } + + /// Get the quadratic phase matrix `B` (cache-aligned). + /// + /// The `±1` phase contributed by `B` is `(-1)^⟨B r, r⟩`. + pub fn aligned_quadratic_phase_matrix(&self) -> &AlignedBitMatrix { + &self.quadratic_phase_matrix + } + + /// Get a copy of the quadratic phase matrix `B` without alignment constraints. + pub fn quadratic_phase_matrix(&self) -> BitMatrix { + BitMatrix::from_aligned(AlignedBitMatrix::from_row_iter( + self.quadratic_phase_matrix.row_iterator(0..self.random_outcome_count()), + self.random_outcome_count(), + )) + } + + /// Get the outcome matrix `M` encoding all `2^{n_r}` measurement branches (cache-aligned). + pub fn aligned_outcome_matrix(&self) -> &AlignedBitMatrix { + &self.outcome_matrix + } + + /// Get a copy of the outcome matrix `M` without alignment constraints. + pub fn outcome_matrix(&self) -> BitMatrix { + BitMatrix::from_aligned(AlignedBitMatrix::from_row_iter( + self.outcome_matrix.row_iterator(0..self.outcome_count()), + self.random_outcome_count(), + )) + } + + /// Get the outcome shift vector `v₀` (cache-aligned). + pub fn aligned_outcome_shift(&self) -> &AlignedBitVec { + &self.outcome_shift + } + + /// Get a copy of the outcome shift vector `v₀` without alignment constraints. + pub fn outcome_shift(&self) -> BitVec { + BitVec::from_aligned(self.outcome_count(), self.outcome_shift.clone()) + } + + /// Get the linear `i`-phase vector `p` (cache-aligned). + /// + /// The `i` phase contributed by `p` is `i^⟨p, r⟩`. + pub fn aligned_linear_i_phase(&self) -> &AlignedBitVec { + &self.linear_i_phase + } + + /// Get a copy of the linear `i`-phase vector `p` without alignment constraints. + pub fn linear_i_phase(&self) -> BitVec { + BitVec::from_aligned(self.random_outcome_count(), self.linear_i_phase.clone()) + } + + /// Get the linear `-1`-phase vector `s` (cache-aligned). + /// + /// The `±1` phase contributed by `s` is `(-1)^⟨s, r⟩`. + pub fn aligned_linear_sign_phase(&self) -> &AlignedBitVec { + &self.linear_sign_phase + } + + /// Get a copy of the linear `-1`-phase vector `s` without alignment constraints. + pub fn linear_sign_phase(&self) -> BitVec { + BitVec::from_aligned(self.random_outcome_count(), self.linear_sign_phase.clone()) + } + + /// Returns the `ζ₈ = e^{iπ/4}` exponent of the scalar `i^⟨p, r⟩ (-1)^⟨B r + s, r⟩` for the given + /// random-bit assignment `r`. + /// + /// This is the `r`-dependent prefactor multiplying `R|A r⟩` in the output state. The remaining + /// phase of `R|A r⟩` itself is carried by [`Self::phased_state_encoder`]. + /// + /// # Panics + /// + /// Panics if `random_bits` has fewer than [`Self::random_outcome_count`] entries. + #[must_use] + pub fn output_phase_exponent(&self, random_bits: &[bool]) -> u8 { + let n_random = self.random_outcome_count(); + assert!( + random_bits.len() >= n_random, + "random_bits is shorter than the number of random outcomes" + ); + phase_form_exponent( + n_random, + |index| random_bits[index], + |index| self.linear_i_phase.index(index), + |index| self.linear_sign_phase.index(index), + |row, column| self.quadratic_phase_matrix.get((row, column)), + ) + } + + /// Sample measurement outcomes from all `2^{n_r}` branches. + pub fn sample(&self, shots: usize) -> BitMatrix { + let mut rng = rand::rng(); + self.sample_with_rng(shots, &mut rng) + } + + /// Sample measurement outcomes using a provided random number generator. + pub fn sample_with_rng(&self, num_shots: usize, rng: &mut R) -> BitMatrix { + let num_outcomes = self.outcome_count(); + let num_random_bits = self.random_outcome_count(); + + if num_outcomes == 0 { + return BitMatrix::from_aligned(AlignedBitMatrix::zeros(num_shots, 0)); + } + + let random_matrix = AlignedBitMatrix::random_with_rng(num_shots, num_random_bits, rng); + let outcome_matrix = + AlignedBitMatrix::from_row_iter(self.outcome_matrix.row_iterator(0..num_outcomes), num_random_bits); + let mut result = random_matrix.mul_transpose(&outcome_matrix); + for shot in 0..num_shots { + result.row_mut(shot).bitxor_assign(&self.outcome_shift.as_view()); + } + BitMatrix::from_aligned(result) + } + + fn ensure_qubit_capacity(&mut self, max_qubit_id: Option) { + if let Some(max_qubit_id) = max_qubit_id { + self.qubit_count = std::cmp::max(self.qubit_count, max_qubit_id + 1); + if max_qubit_id >= self.qubit_capacity() { + let new_capacity = (max_qubit_id + 1).next_power_of_two(); + self.reserve_qubits(new_capacity); + } + } + } + + #[inline] + fn ensure_outcome_capacity(&mut self, random_outcome: bool) { + let mut new_outcome_capacity = self.outcome_capacity(); + let next_outcome_pos = self.outcome_count(); + if next_outcome_pos >= self.outcome_capacity() { + new_outcome_capacity = (next_outcome_pos + 1).next_power_of_two(); + } + + let mut new_random_outcome_capacity = self.random_outcome_capacity(); + if random_outcome { + let next_random_bit = self.random_outcome_count(); + if next_random_bit >= self.random_outcome_capacity() { + new_random_outcome_capacity = (next_random_bit + 1).next_power_of_two(); + } + } + + self.reserve_outcomes(new_outcome_capacity, new_random_outcome_capacity); + } + + /// Applies a Pauli `P` conditioned on the random-bit parity `⟨indicator, r⟩`, updating the + /// state matrix `A` and the phase data `B`, `p`, `s` accordingly. + /// + /// This realises step 4b of Algorithm 4.2: with preimage `i^l X^x Z^z = R† P R`, + /// `A ← A + x indicator^T` and the conditional factor `(i^l (-1)^{⟨z, A r⟩})^⟨indicator, r⟩` is + /// absorbed into the tracked phase: + /// + /// * the `(-1)^{⟨z, A r⟩ ⟨indicator, r⟩}` factor adds `indicator (z^T A)^T` to `B`; + /// * for odd `l`, the extra `i^⟨indicator, r⟩` flips `p` by `indicator`, and its carry against the + /// existing i-phase, `(-1)^{⟨p, r⟩⟨indicator, r⟩}`, adds the outer product `p ⊗ indicator` to `B`; + /// * for `l` with the two-bit set, the `(-1)^⟨indicator, r⟩` factor flips `s` by `indicator`. + fn apply_pauli_conditioned_on_inner_random_bits( + &mut self, + pauli: &PauliUnitary, + inner_bits_indicator: &AlignedBitVec, + ) { + let preimage = self.phased_clifford.clifford().preimage(pauli); + let z_times_sign_matrix = row_sum(&self.sign_matrix, preimage.z_bits().support()); + for row in inner_bits_indicator.support() { + self.quadratic_phase_matrix + .row_mut(row) + .bitxor_assign(&z_times_sign_matrix); + } + for x_bit_pos in preimage.x_bits().support() { + self.sign_matrix.row_mut(x_bit_pos).bitxor_assign(inner_bits_indicator); + } + let l = preimage.xz_phase_exponent().value(); + if l & 1 == 1 { + // Adding i^⟨indicator, r⟩ to the i-phase, which is tracked mod 2 in `p`. The carry + // i·i = -1 from overlap with the existing i-phase is the quadratic cross term + // (-1)^{⟨p, r⟩⟨indicator, r⟩}, absorbed into `B` as the outer product `p ⊗ indicator`. + for row in self.linear_i_phase.support() { + self.quadratic_phase_matrix + .row_mut(row) + .bitxor_assign(inner_bits_indicator); + } + self.linear_i_phase.bitxor_assign(inner_bits_indicator); + } + if l & 2 == 2 { + self.linear_sign_phase.bitxor_assign(inner_bits_indicator); + } + } + + pub fn with_capacity(qubit_count: usize, outcome_count: usize, random_outcome_count: usize) -> Self { + const MIN_CAPACITY: usize = CLIFFORD_BIT_ALIGNMENT; + let outcome_capacity = outcome_count.max(MIN_CAPACITY); + let random_capacity = random_outcome_count.max(MIN_CAPACITY); + + PhasedOutcomeCompleteSimulation { + phased_clifford: PhasedCliffordUnitary::identity(qubit_count), + sign_matrix: AlignedBitMatrix::zeros(qubit_count, random_capacity), + quadratic_phase_matrix: AlignedBitMatrix::zeros(random_capacity, random_capacity), + outcome_matrix: AlignedBitMatrix::zeros(outcome_capacity, random_capacity), + outcome_shift: AlignedBitVec::zeros(outcome_capacity), + linear_i_phase: AlignedBitVec::zeros(random_capacity), + linear_sign_phase: AlignedBitVec::zeros(random_capacity), + random_outcome_indicator: Vec::with_capacity(outcome_count), + symbolic_angle_indicator: Vec::with_capacity(random_outcome_count), + random_bit_count: 0, + qubit_count, + } + } + + /// Measures a Pauli observable using an anti-commuting hint operator, tracking the exact phase. + /// + /// Implements case 5 of Algorithm 4.2. Given an anti-commuting hint `P'` with preimage + /// `R† P' R = (-1)^α Z^{b'}`, the encoder is updated by `R ← e^{iπ/4 (i P' P)} R`, the quadratic + /// and linear `-1` phases absorb the outcome-dependent stabiliser sign, and the `(-1)^α` sign + /// relabels the reported outcome (`m = r ⊕ α`) via `outcome_shift` rather than a global phase. + /// + /// # Panics + /// + /// Panics if `hint` does not anti-commute with `observable`. + pub fn measure_pauli_with_hint_generic( + &mut self, + observable: &SparsePauli, + hint: &PauliUnitary, + ) { + assert!( + anti_commutes_with(observable, hint), + "observable={observable}, hint={hint}" + ); + let preimage = self.phased_clifford.clifford().preimage(hint); + if preimage.x_bits().support().next().is_some() { + // hint is not a stabilizer of the encoded state family + self.measure(observable); + } else { + // Ensure capacity for the new random bit before sizing the indicator vectors. + self.ensure_outcome_capacity(true); + + // R <- (-1)^alpha e^{i pi/4 (i P' P)} R. + // i P' P = -i P P', so the rotation Pauli is (observable * hint) with an i^3 = -i phase. + let alpha = preimage.xz_phase_exponent().value() / 2; + let mut rotation = observable.clone(); + rotation.mul_assign_right(hint); + rotation.add_assign_phase_exp(3); + self.phased_clifford.left_mul_pauli_exp(&rotation); + + // a = A^T b', with the new random bit appended: a_with_zero and a_with_one = a ⊕ {0,1}. + let a_with_zero = row_sum(&self.sign_matrix, preimage.z_bits().support()); + let mut a_with_one = a_with_zero.clone(); + a_with_one.assign_index(self.random_bit_count, true); + let new_random_bit = self.random_bit_count; + self.allocate_random_bit(); + + // B <- B + (a ⊕ 0)(a ⊕ 1)^T, s_{n(s)} <- s_{n(s)} + alpha. + for row in a_with_zero.support() { + self.quadratic_phase_matrix.row_mut(row).bitxor_assign(&a_with_one); + } + if alpha == 1 { + self.linear_sign_phase + .assign_index(new_random_bit, !self.linear_sign_phase.index(new_random_bit)); + // (-1)^alpha relabels the reported outcome (m = r ⊕ alpha), not the global phase. + let outcome_position = self.outcome_count() - 1; + self.outcome_shift + .assign_index(outcome_position, !self.outcome_shift.index(outcome_position)); + } + + // Apply P' conditioned on the random bits indicated by (a ⊕ 1). + self.apply_pauli_conditioned_on_inner_random_bits(hint, &a_with_one); + } + } + + fn measure_deterministic(&mut self, preimage: &PauliUnitary) { + self.ensure_outcome_capacity(false); + let outcome_matrix_row = row_sum(&self.sign_matrix, preimage.z_bits().support()); + let outcome_position = self.random_outcome_indicator.len(); + self.outcome_matrix + .row_mut(outcome_position) + .assign(&outcome_matrix_row); + debug_assert!(preimage.xz_phase_exponent().is_even()); + if preimage.xz_phase_exponent().value() == 2 { + self.outcome_shift.assign_index(outcome_position, true); + } + self.random_outcome_indicator.push(false); + } + + /// Get the number of random (non-deterministic) measurement outcomes. + #[must_use] + pub fn random_outcome_count(&self) -> usize { + self.random_bit_count + } + + /// Get indicators for which outcomes are random. + #[must_use] + pub fn random_outcome_indicator(&self) -> &[bool] { + &self.random_outcome_indicator + } + + /// Get indicators for which random bits are symbolic rotation angles. + /// + /// The returned slice is indexed by random-bit index (`0..random_outcome_count()`). Entry `k` is + /// `true` iff random bit `k` was allocated via [`Simulation::allocate_symbolic_angle`] (a virtual + /// rotation parameter) rather than [`Simulation::allocate_random_bit`] or a genuine measurement. + #[must_use] + pub fn symbolic_angle_indicator(&self) -> &[bool] { + &self.symbolic_angle_indicator + } + + fn allocate_random_bit_with_provenance(&mut self, is_symbolic_angle: bool) -> usize { + self.ensure_outcome_capacity(true); + let outcome_pos = self.random_outcome_indicator.len(); + self.outcome_matrix + .row_mut(outcome_pos) + .assign_index(self.random_bit_count, true); + self.random_outcome_indicator.push(true); + self.symbolic_angle_indicator.push(is_symbolic_angle); + self.random_bit_count += 1; + self.random_bit_count - 1 + } +} + +impl Simulation for PhasedOutcomeCompleteSimulation { + fn allocate_random_bit(&mut self) -> usize { + self.allocate_random_bit_with_provenance(false) + } + + fn allocate_symbolic_angle(&mut self) -> usize { + self.allocate_random_bit_with_provenance(true) + } + + fn clifford(&mut self, _clifford: &crate::Unitary, _support: &[crate::QubitId]) { + unimplemented!( + "PhasedOutcomeCompleteSimulation tracks the exact global phase, which a phaseless \ + CliffordUnitary does not determine; apply Cliffords through unitary_op, pauli or \ + pauli_exp instead" + ); + } + + fn unitary_op(&mut self, unitary_op: UnitaryOp, support: &[crate::QubitId]) { + self.ensure_qubit_capacity(max_support(support)); + self.phased_clifford.left_mul(unitary_op, support); + } + + fn permute(&mut self, permutation: &[usize], support: &[crate::QubitId]) { + self.ensure_qubit_capacity(max_support(support)); + self.phased_clifford.left_mul_permutation(permutation, support); + } + + fn controlled_pauli(&mut self, observable1: &SparsePauli, observable2: &SparsePauli) { + self.ensure_qubit_capacity(max_pair_support(observable1, observable2)); + self.controlled_pauli_phase_resolved(observable1, observable2); + } + + fn pauli(&mut self, observable: &SparsePauli) { + self.ensure_qubit_capacity(observable.max_support()); + self.phased_clifford.left_mul_pauli(observable); + } + + fn pauli_exp(&mut self, observable: &SparsePauli) { + self.ensure_qubit_capacity(observable.max_support()); + self.phased_clifford.left_mul_pauli_exp(observable); + } + + fn is_stabilizer_up_to_sign(&self, observable: &SparsePauli) -> bool { + self.phased_clifford.clifford().preimage(observable).x_bits().is_zero() + } + + fn qubit_count(&self) -> usize { + self.qubit_count + } + + fn conditional_pauli(&mut self, observable: &SparsePauli, outcomes: &[usize], parity: bool) { + self.ensure_qubit_capacity(observable.max_support()); + let bit_indicator = outcomes.iter().copied().collect::(); + let is_p_applied: bool = !parity ^ bit_indicator.dot(&self.outcome_shift); + if is_p_applied { + self.pauli(observable); + } + let inner_bits_indicator = row_sum(&self.outcome_matrix, outcomes); + self.apply_pauli_conditioned_on_inner_random_bits(observable, &inner_bits_indicator); + } + + fn is_stabilizer(&self, observable: &SparsePauli) -> bool { + let preimage = self.phased_clifford.clifford().preimage(observable); + if preimage.x_bits().is_zero() { + let sign_parity_indicator = row_sum(&self.sign_matrix, preimage.z_bits().support()); + sign_parity_indicator.is_zero() + } else { + false + } + } + + fn is_stabilizer_with_conditional_sign(&self, observable: &SparsePauli, outcomes: &[crate::OutcomeId]) -> bool { + let preimage = self.phased_clifford.clifford().preimage(observable); + if preimage.x_bits().is_zero() { + let sign_parity_indicator = row_sum(&self.sign_matrix, preimage.z_bits().support()); + debug_assert!(preimage.xz_phase_exponent().is_even()); + let shift = preimage.xz_phase_exponent().value() / 2 == 1; + let expected_parity_indicator = row_sum(&self.outcome_matrix, outcomes.iter().copied()); + let expected_shift = outcomes + .iter() + .copied() + .map(|o| self.outcome_shift.index(o)) + .fold(false, |acc, v| acc ^ v); + (sign_parity_indicator == expected_parity_indicator) && (shift == expected_shift) + } else { + false + } + } + + fn measure(&mut self, observable: &SparsePauli) -> usize { + self.ensure_qubit_capacity(observable.max_support()); + let preimage = self.phased_clifford.clifford().preimage(observable); + let non_zero_pos = preimage.x_bits().support().next(); + match non_zero_pos { + Some(pos) => { + let hint = self.phased_clifford.clifford().image_z(pos); + self.measure_pauli_with_hint_generic(observable, &hint); + } + None => { + self.measure_deterministic(&preimage); + } + } + self.outcome_count() - 1 + } + + fn measure_with_hint(&mut self, observable: &SparsePauli, hint: &SparsePauli) -> usize { + self.ensure_qubit_capacity(max_pair_support(observable, hint)); + self.measure_pauli_with_hint_generic(observable, hint); + self.outcome_count() - 1 + } + + fn outcome_count(&self) -> usize { + self.random_outcome_indicator.len() + } + + fn with_capacity(qubit_count: usize, outcome_count: usize, random_outcome_count: usize) -> Self + where + Self: Sized, + { + PhasedOutcomeCompleteSimulation::with_capacity(qubit_count, outcome_count, random_outcome_count) + } + + fn qubit_capacity(&self) -> usize { + debug_assert_eq!(self.phased_clifford.num_qubits(), self.sign_matrix.row_count()); + self.phased_clifford.num_qubits() + } + + fn outcome_capacity(&self) -> usize { + self.outcome_matrix.row_count() + } + + fn random_outcome_capacity(&self) -> usize { + debug_assert_eq!(self.outcome_matrix.column_count(), self.sign_matrix.column_count()); + self.outcome_matrix.column_count() + } + + fn reserve_qubits(&mut self, new_capacity: usize) { + if new_capacity > self.qubit_capacity() { + self.sign_matrix.resize(new_capacity, self.sign_matrix.column_count()); + self.phased_clifford.resize(new_capacity); + } + } + + fn reserve_outcomes(&mut self, new_outcome_capacity: usize, new_random_outcome_capacity: usize) { + assert!( + new_outcome_capacity >= new_random_outcome_capacity, + "outcome capacity must be at least random outcome capacity" + ); + let new_outcome_capacity = new_outcome_capacity.max(self.outcome_capacity()); + let new_random_outcome_capacity = new_random_outcome_capacity.max(self.random_outcome_capacity()); + + self.outcome_matrix + .resize(new_outcome_capacity, new_random_outcome_capacity); + self.sign_matrix + .resize(self.sign_matrix.row_count(), new_random_outcome_capacity); + self.quadratic_phase_matrix + .resize(new_random_outcome_capacity, new_random_outcome_capacity); + if self.outcome_shift.len() < new_outcome_capacity { + self.outcome_shift.resize(new_outcome_capacity); + } + if self.linear_i_phase.len() < new_random_outcome_capacity { + self.linear_i_phase.resize(new_random_outcome_capacity); + } + if self.linear_sign_phase.len() < new_random_outcome_capacity { + self.linear_sign_phase.resize(new_random_outcome_capacity); + } + } +} + +impl PhasedOutcomeCompleteSimulation { + /// Applies the controlled-Pauli `Λ(observable1, observable2)` to the encoder, phase-resolved. + /// + /// For commuting Hermitian involutions `P₁`, `P₂` the controlled-Pauli is + /// `Λ(P₁, P₂) = (I + P₁)/2 + (I − P₁)/2 · P₂ = exp(iπ/4 · (I − P₁)(I − P₂))`. Since the factors + /// commute this is + /// + /// ```text + /// Λ(P₁, P₂) = e^{iπ/4} · e^{-iπ/4 P₁} · e^{-iπ/4 P₂} · e^{iπ/4 P₁P₂}, + /// ``` + /// + /// each factor of which is applied through the phase-exact primitive, so the global phase of the + /// encoder is tracked exactly while the symplectic action matches + /// [`paulimer::clifford::CliffordMutable::left_mul_controlled_pauli`]. + fn controlled_pauli_phase_resolved(&mut self, observable1: &SparsePauli, observable2: &SparsePauli) { + debug_assert!( + paulimer::pauli::commutes_with(observable1, observable2), + "controlled_pauli requires commuting observables" + ); + let mut negated1 = observable1.clone(); + negated1.add_assign_phase_exp(2); + let mut negated2 = observable2.clone(); + negated2.add_assign_phase_exp(2); + let mut product = observable1.clone(); + product.mul_assign_right(observable2); + + self.phased_clifford.left_mul_pauli_exp(&negated1); + self.phased_clifford.left_mul_pauli_exp(&negated2); + self.phased_clifford.left_mul_pauli_exp(&product); + self.phased_clifford.left_mul_global_phase(1); + } +} diff --git a/pauliverse/tests/action_test.rs b/pauliverse/tests/action_test.rs index 44cd6f27..8b250832 100644 --- a/pauliverse/tests/action_test.rs +++ b/pauliverse/tests/action_test.rs @@ -98,6 +98,13 @@ proptest! { check_and_compare_unitary(&z_diagonal_unitary, &circuit, &input); } + #[test] + fn diagonal_unitary_x_ejection_proptest(z_diagonal_unitary in arbitrary_diagonal_clifford(1..6usize)) { + let x_diagonal_unitary = x_diagonal_from_z_diagonal(&z_diagonal_unitary); + let (circuit, input) = diagonal_unitary_x_ejection_circuit_with_io(&x_diagonal_unitary); + check_and_compare_unitary_x(&x_diagonal_unitary, &circuit, &input); + } + #[test] fn diagonal_unitary_injection_proptest(z_diagonal_unitary in arbitrary_diagonal_clifford(1..6usize)) { let (circuit, input) = diagonal_unitary_injection_circuit_with_io(&z_diagonal_unitary); @@ -136,6 +143,37 @@ proptest! { ); } + #[test] + fn diagonal_measure_x_ejection_proptest(z_diagonal_paulis in arbitrary_independent_z_paulis(1..6usize, 1..3usize)) { + let x_diagonal_paulis = x_paulis_from_z_paulis(&z_diagonal_paulis); + let (ejection_circuit, input_output, outcome_map) = diagonal_measure_x_ejection_circuit_with_io(&x_diagonal_paulis); + let ejection_action = + action_of(&ejection_circuit, &input_output, &input_output).expect("diagonal X measure ejection action"); + + check_multi_pauli_action(&x_diagonal_paulis, &input_output, &input_output, &ejection_action); + + let (measure_circuit, measure_input_output) = multi_measure_circuit_with_io(&x_diagonal_paulis); + let measure_action = + action_of(&measure_circuit, &measure_input_output, &measure_input_output).expect("diagonal X measure action"); + + check_multi_pauli_action( + &x_diagonal_paulis, + &measure_input_output, + &measure_input_output, + &measure_action, + ); + let map = affine_map_from_sparse( + ejection_action.outcome_count(), + measure_action.outcome_count(), + outcome_map, + ); + measure_action + .is_equivalent_with_map(&ejection_action, Some(&map)) + .expect( + "diagonal X measure ejection action should be equivalent to diagonal X measure action with outcome mapping", + ); + } + #[test] fn diagonal_measure_injection_proptest(z_diagonal_paulis in arbitrary_independent_z_paulis(1..6usize, 1..3usize)) { let (injection_circuit, input_output, _) = diagonal_measure_injection_circuit_with_io(&z_diagonal_paulis); @@ -177,6 +215,29 @@ fn check_and_compare_unitary( .expect("diagonal ejection action should be equivalent to unitary action"); } +fn check_and_compare_unitary_x( + x_diagonal_unitary: &CliffordUnitary, + circuit: &Circuit, + input_and_output_qubits: &[usize], +) { + assert!(x_diagonal_unitary.is_diagonal(XOrZ::X)); + let action = + action_of(circuit, input_and_output_qubits, input_and_output_qubits).expect("X diagonal ejection action"); + check_unitary_action( + x_diagonal_unitary, + input_and_output_qubits, + input_and_output_qubits, + &action, + ); + + let (unitary_circuit, unitary_input, unitary_output) = one_unitary_circuit_with_io(x_diagonal_unitary); + let unitary_action = + action_of(&unitary_circuit, &unitary_input, &unitary_output).expect("X diagonal unitary action"); + unitary_action + .is_equivalent_with_map(&action, None) + .expect("X diagonal ejection action should be equivalent to unitary action"); +} + /// Validation of some common kinds of action fn check_bell_pair(circuit: &Circuit, input: &[usize], output: &[usize]) { let action = action_of(circuit, input, output).expect("Bell pair preparation action"); @@ -444,6 +505,50 @@ fn diagonal_unitary_ejection_circuit_with_io(z_diagonal_unitary: &CliffordUnitar (b.into_circuit(), targets) } +/// A transversal Hadamard layer on `qubit_count` qubits, used to map between the Z- and X-diagonal +/// Clifford subgroups by conjugation. +fn hadamard_layer(qubit_count: usize) -> CliffordUnitary { + let mut layer = CliffordUnitary::identity(qubit_count); + for qubit in 0..qubit_count { + layer.left_mul(UnitaryOp::Hadamard, &[qubit]); + } + layer +} + +/// Conjugates a Z-diagonal Clifford by a transversal Hadamard to obtain an X-diagonal Clifford +/// `H^n · U · H^n`. (The result is independent of association since `H` is its own inverse.) +fn x_diagonal_from_z_diagonal(z_diagonal_unitary: &CliffordUnitary) -> CliffordUnitary { + let hadamards = hadamard_layer(z_diagonal_unitary.num_qubits()); + &(&hadamards * z_diagonal_unitary) * &hadamards +} + +/// Implements `x_diagonal_unitary` via the X-basis dual of the diagonal ejection of Figure 9 in +/// : the whole gadget is the conjugation of +/// [`diagonal_unitary_ejection_circuit_with_io`] by a transversal Hadamard on every system and +/// reference qubit. Each reference (ancilla) is prepared in `|+⟩`, the CNOTs run from references into +/// the targets, an X-diagonal Clifford acts on the references, and the references are measured in the +/// Z basis with a conditional X correction on a `1` outcome. +fn diagonal_unitary_x_ejection_circuit_with_io(x_diagonal_unitary: &CliffordUnitary) -> (Circuit, Vec) { + assert!(x_diagonal_unitary.is_diagonal(XOrZ::X)); + let qubit_count = x_diagonal_unitary.num_qubits(); + let targets = (0..qubit_count).collect::>(); + let references = (qubit_count..2 * qubit_count).collect::>(); + + let mut b = empty_builder(); + for &reference in &references { + b = b.h(reference); + } + for (&target, &reference) in targets.iter().zip(references.iter()) { + b = b.cnot(reference, target); + } + b = b.clifford(x_diagonal_unitary, &references); + for (id, (&target, &reference)) in targets.iter().zip(references.iter()).enumerate() { + b = b.measure_z(reference, id).conditional_x(target, &[id], true); + } + + (b.into_circuit(), targets) +} + type OutcomeMapping = Vec<(OutcomeId, bool, Vec)>; /// Implements measurement of `z_diagonal_paulis` via diagonal ejection, see Figure 9 in @@ -477,6 +582,54 @@ fn diagonal_measure_ejection_circuit_with_io( (b.into_circuit(), targets, outcome_map) } +/// X-basis dual of [`x_paulis_from_z_paulis`]'s input: turns each Z-Pauli into the X-Pauli with the +/// same support (swapping the `Z` and `X` bits), so independent Z-Paulis become independent X-Paulis. +fn x_paulis_from_z_paulis(z_diagonal_paulis: &[SparsePauli]) -> Vec { + z_diagonal_paulis + .iter() + .map(|pauli| SparsePauli::from_bits(pauli.z_bits().clone(), IndexSet::new(), 0)) + .collect() +} + +/// Implements measurement of `x_diagonal_paulis` via the X-basis dual of the diagonal measure +/// ejection of Figure 9 in : the whole gadget is the conjugation +/// of [`diagonal_measure_ejection_circuit_with_io`] by a transversal Hadamard. Each reference is +/// prepared in `|+⟩`, the CNOTs run from references into the targets, the X-diagonal Paulis are +/// measured (non-destructively) on the references, and each reference is destructively measured in +/// the Z basis with a conditional X correction on a `1` outcome. +fn diagonal_measure_x_ejection_circuit_with_io( + x_diagonal_paulis: &[SparsePauli], +) -> (Circuit, Vec, OutcomeMapping) { + let qubit_count = max_qubit_id_of(x_diagonal_paulis) + 1; + + let targets = (0..qubit_count).collect::>(); + let references = (qubit_count..2 * qubit_count).collect::>(); + + let mut b = empty_builder(); + for &reference in &references { + b = b.h(reference); + } + for (&target, &reference) in targets.iter().zip(references.iter()) { + b = b.cnot(reference, target); + } + + let next_outcome_id = b.next_outcome_id(); + let pauli_outcome_ids = next_outcome_id..(next_outcome_id + x_diagonal_paulis.len()); + for (pauli, outcome_id) in x_diagonal_paulis.iter().zip(pauli_outcome_ids.clone()) { + let reference_pauli = remapped_sparse(pauli, &references); + b = b.measure_sparse(&reference_pauli, outcome_id); + } + + let next_outcome_id = b.next_outcome_id(); + let z_outcome_ids = next_outcome_id..(next_outcome_id + targets.len()); + for (id, (&target, &reference)) in z_outcome_ids.zip(targets.iter().zip(references.iter())) { + b = b.measure_z(reference, id).conditional_x(target, &[id], true); + } + + let outcome_map = pauli_outcome_ids.map(|id| (id, false, vec![id])).collect::>(); + (b.into_circuit(), targets, outcome_map) +} + fn multi_measure_circuit_with_io(z_diagonal_paulis: &[SparsePauli]) -> (Circuit, Vec) { let max_qubit_id = max_qubit_id_of(z_diagonal_paulis); let qubits = (0..=max_qubit_id).collect::>(); diff --git a/pauliverse/tests/measure_with_hint_sign_test.rs b/pauliverse/tests/measure_with_hint_sign_test.rs new file mode 100644 index 00000000..71383985 --- /dev/null +++ b/pauliverse/tests/measure_with_hint_sign_test.rs @@ -0,0 +1,94 @@ +//! Sign-correctness of hinted Pauli measurements (`measure_with_hint`). +//! +//! Defining property of a projective measurement: once `P` has been measured, `P` is a stabilizer of +//! the post-measurement state whose sign equals the reported outcome. In particular, this must hold +//! no matter which anti-commuting `hint` is supplied — including hints carrying a negative sign, +//! which drive the `(-1)^alpha` branch of case 5 of Algorithm 4.2. + +use paulimer::UnitaryOp; +use paulimer::clifford::{Clifford, PhasedCliffordUnitary}; +use paulimer::pauli::{SparsePauli, as_sparse}; +use pauliverse::{PhasedOutcomeCompleteSimulation, Simulation}; +use proptest::collection::vec; +use proptest::prelude::*; + +#[derive(Clone, Debug)] +enum Gate { + Single { op: UnitaryOp, qubit: usize }, + Two { op: UnitaryOp, first: usize, second: usize }, +} + +fn distinct_pair(qubit_count: usize) -> impl Strategy { + (0..qubit_count, 0..qubit_count - 1) + .prop_map(|(first, second)| (first, if second < first { second } else { second + 1 })) +} + +fn gate_strategy(qubit_count: usize) -> BoxedStrategy { + use UnitaryOp::{ControlledX, Hadamard, SqrtX, SqrtZ}; + let single = (prop::sample::select(vec![Hadamard, SqrtZ, SqrtX]), 0..qubit_count) + .prop_map(|(op, qubit)| Gate::Single { op, qubit }); + if qubit_count >= 2 { + let two = distinct_pair(qubit_count).prop_map(|(first, second)| Gate::Two { + op: ControlledX, + first, + second, + }); + prop_oneof![3 => single, 1 => two].boxed() + } else { + single.boxed() + } +} + +/// Generates a qubit count, a random Clifford preparation circuit, whether to negate the hint, and +/// the qubit whose stabilizer/destabilizer images drive the measurement. +fn scenario() -> impl Strategy, bool, usize)> { + (1usize..4).prop_flat_map(|qubit_count| { + ( + Just(qubit_count), + vec(gate_strategy(qubit_count), 0..12), + any::(), + 0..qubit_count, + ) + }) +} + +fn apply(gate: &Gate, sim: &mut PhasedOutcomeCompleteSimulation, mirror: &mut PhasedCliffordUnitary) { + match *gate { + Gate::Single { op, qubit } => { + sim.unitary_op(op, &[qubit]); + mirror.left_mul(op, &[qubit]); + } + Gate::Two { op, first, second } => { + sim.unitary_op(op, &[first, second]); + mirror.left_mul(op, &[first, second]); + } + } +} + +proptest! { + /// For a random stabilizer state, measuring the destabilizer `X`-image of qubit `target` while + /// hinting with the (optionally negated) stabilizer `Z`-image of `target` must leave the + /// observable a stabilizer whose conditional sign matches the reported outcome. + #[test] + fn measure_with_hint_outcome_sign_is_correct((qubit_count, gates, negate_hint, target) in scenario()) { + let mut sim = PhasedOutcomeCompleteSimulation::new(qubit_count); + let mut mirror = PhasedCliffordUnitary::identity(qubit_count); + for gate in &gates { + apply(gate, &mut sim, &mut mirror); + } + + // `image_z(target)` is a stabilizer of the prepared state; `image_x(target)` anti-commutes + // with it, so measuring the latter is a genuine (random) case-5 measurement. + let stabilizer: SparsePauli = as_sparse(&mirror.clifford().image_z(target)); + let observable: SparsePauli = as_sparse(&mirror.clifford().image_x(target)); + let hint = if negate_hint { -stabilizer } else { stabilizer }; + + let outcome = sim.measure_with_hint(&observable, &hint); + + prop_assert!( + sim.is_stabilizer_with_conditional_sign(&observable, &[outcome]), + "measured observable {observable} is not a stabilizer with the reported outcome sign \ + (qubit_count={qubit_count}, negate_hint={negate_hint}, target={target})" + ); + } +} diff --git a/pauliverse/tests/phased_action_test.rs b/pauliverse/tests/phased_action_test.rs new file mode 100644 index 00000000..ae4f8449 --- /dev/null +++ b/pauliverse/tests/phased_action_test.rs @@ -0,0 +1,959 @@ +use paulimer::core::{x, z}; +use paulimer::pauli::SparsePauli; +use paulimer::{PositionedPauliObservable, UnitaryOp}; +use pauliverse::action::{ + ActionsInequivalenceReason, PhasedCircuitAction, phased_action_from_simulation, phased_action_of, +}; +use pauliverse::phased_outcome_complete_simulation::PhasedOutcomeCompleteSimulation; +use pauliverse::{Circuit, CircuitBuilder, QubitId, Simulation}; +use proptest::prelude::*; +use rand::SeedableRng; +use std::ops::Range; + +fn build_circuit(build: impl FnOnce(&mut CircuitBuilder)) -> Circuit { + let mut builder = CircuitBuilder::new(); + build(&mut builder); + builder.into() +} + +fn sparse(observable: &[PositionedPauliObservable]) -> SparsePauli { + observable.into() +} + +/// `exp(iα Z₀Z₁)` represented as a symbolic rotation gadget: allocate a random branch bit, then +/// conditionally apply `Z₀Z₁` on the odd branch. +fn zz_rotation() -> (Circuit, Vec, Vec) { + let circuit = build_circuit(|builder| { + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0), z(1)]), branch); + }); + (circuit, vec![0, 1], vec![0, 1]) +} + +/// `CNOT₀₁ · exp(iα Z₁) · CNOT₀₁`, which should equal `exp(iα Z₀Z₁)` as a channel. +fn cnot_conjugated_z_rotation() -> (Circuit, Vec, Vec) { + let circuit = build_circuit(|builder| { + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(1)]), branch); + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + }); + (circuit, vec![0, 1], vec![0, 1]) +} + +/// `exp(iα Z₁)` on its own, which differs from `exp(iα Z₀Z₁)` in symplectic action. +fn z_rotation() -> (Circuit, Vec, Vec) { + let circuit = build_circuit(|builder| { + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(1)]), branch); + }); + (circuit, vec![0, 1], vec![0, 1]) +} + +/// Conditional `±Z₀` gadget: the two signs share the same symplectic action but differ only in the +/// branch phase of the odd branch. +fn signed_z_rotation(negate: bool) -> (Circuit, Vec, Vec) { + let circuit = build_circuit(|builder| { + let branch = builder.allocate_symbolic_angle(); + let observable = if negate { -sparse(&[z(0)]) } else { sparse(&[z(0)]) }; + builder.symbolic_pauli_exp(&observable, branch); + }); + (circuit, vec![0], vec![0]) +} + +#[test] +fn zz_rotation_equals_cnot_conjugated_z_rotation() { + let (direct, direct_input, direct_output) = zz_rotation(); + let (conjugated, conjugated_input, conjugated_output) = cnot_conjugated_z_rotation(); + + let direct_action = phased_action_of(&direct, &direct_input, &direct_output).expect("direct action"); + let conjugated_action = + phased_action_of(&conjugated, &conjugated_input, &conjugated_output).expect("conjugated action"); + + direct_action + .is_equivalent(&conjugated_action) + .expect("symbolic rotations must agree including branch phase"); + conjugated_action + .is_equivalent(&direct_action) + .expect("equivalence must be symmetric"); +} + +#[test] +fn zz_rotation_differs_from_z_rotation() { + let (zz, zz_input, zz_output) = zz_rotation(); + let (single, single_input, single_output) = z_rotation(); + + let zz_action = phased_action_of(&zz, &zz_input, &zz_output).expect("zz action"); + let single_action = phased_action_of(&single, &single_input, &single_output).expect("z action"); + + let reasons = zz_action + .is_equivalent(&single_action) + .expect_err("rotations with different supports must differ"); + assert!(!reasons.is_empty()); +} + +#[test] +fn opposite_sign_rotations_differ_only_in_relative_phase() { + let (positive, positive_input, positive_output) = signed_z_rotation(false); + let (negative, negative_input, negative_output) = signed_z_rotation(true); + + let positive_action = phased_action_of(&positive, &positive_input, &positive_output).expect("positive action"); + let negative_action = phased_action_of(&negative, &negative_input, &negative_output).expect("negative action"); + + positive_action + .is_equivalent_up_to_signs(&negative_action) + .expect("phaseless actions must be identical"); + + let reasons = positive_action + .is_equivalent(&negative_action) + .expect_err("opposite signs must be distinguished by the phased action"); + assert_eq!(reasons, vec![ActionsInequivalenceReason::RelativePhase]); +} + +#[test] +fn rotation_equals_itself() { + let (zz, zz_input, zz_output) = zz_rotation(); + let action = phased_action_of(&zz, &zz_input, &zz_output).expect("action"); + action + .is_equivalent(&action) + .expect("a rotation must be equivalent to itself"); +} + +/// Conjugating a `Z` rotation by a Hadamard produces the corresponding `X` rotation: +/// `H₀ · exp(iα Z₀) · H₀ == exp(iα X₀)`, including branch phase. +#[test] +fn hadamard_conjugated_z_rotation_equals_x_rotation() { + let conjugated = build_circuit(|builder| { + builder.unitary_op(UnitaryOp::Hadamard, &[0]); + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0)]), branch); + builder.unitary_op(UnitaryOp::Hadamard, &[0]); + }); + let x_rotation = build_circuit(|builder| { + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[x(0)]), branch); + }); + + let conjugated_action = phased_action_of(&conjugated, &[0], &[0]).expect("conjugated action"); + let x_action = phased_action_of(&x_rotation, &[0], &[0]).expect("x action"); + + conjugated_action + .is_equivalent(&x_action) + .expect("HZH must equal an X rotation, including branch phase"); + x_action + .is_equivalent(&conjugated_action) + .expect("equivalence must be symmetric"); +} + +/// A mixed-basis exponent equals its single-qubit Hadamard conjugate: +/// `exp(iα X₀Z₁) == H₀ · exp(iα Z₀Z₁) · H₀`, while the un-conjugated `exp(iα Z₀Z₁)` differs. +#[test] +fn mixed_basis_exponent_equals_hadamard_conjugated_zz() { + let mixed = build_circuit(|builder| { + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[x(0), z(1)]), branch); + }); + let conjugated_zz = build_circuit(|builder| { + builder.unitary_op(UnitaryOp::Hadamard, &[0]); + let branch = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0), z(1)]), branch); + builder.unitary_op(UnitaryOp::Hadamard, &[0]); + }); + + let mixed_action = phased_action_of(&mixed, &[0, 1], &[0, 1]).expect("mixed action"); + let conjugated_action = phased_action_of(&conjugated_zz, &[0, 1], &[0, 1]).expect("conjugated action"); + + mixed_action + .is_equivalent(&conjugated_action) + .expect("mixed-basis exponent must equal the Hadamard-conjugated ZZ exponent"); + + let (bare_zz, bare_input, bare_output) = zz_rotation(); + let bare_action = phased_action_of(&bare_zz, &bare_input, &bare_output).expect("bare zz action"); + mixed_action + .is_equivalent(&bare_action) + .expect_err("the un-conjugated ZZ exponent is a different channel"); +} + +/// Builds the Choi state of a single-system-qubit gadget directly in a phased simulation, mirroring +/// the simulator-native idiom used by the Python bindings: Bell-pair every system qubit `q` in +/// `0..n` with its reference `q + n`, allocate one random branch bit, then apply `build_gadget`. +fn choi_simulation( + system_qubit_count: usize, + build_gadget: impl FnOnce(&mut PhasedOutcomeCompleteSimulation, usize), +) -> PhasedOutcomeCompleteSimulation { + let mut simulation = PhasedOutcomeCompleteSimulation::new(2 * system_qubit_count); + for system_qubit in 0..system_qubit_count { + simulation.unitary_op( + UnitaryOp::PrepareBell, + &[system_qubit, system_qubit + system_qubit_count], + ); + } + let branch = simulation.allocate_symbolic_angle(); + build_gadget(&mut simulation, branch); + simulation +} + +#[test] +fn simulator_native_action_matches_circuit_action() { + let (zz, zz_input, zz_output) = zz_rotation(); + let circuit_action = phased_action_of(&zz, &zz_input, &zz_output).expect("circuit action"); + + let simulation = choi_simulation(2, |simulation, branch| { + simulation.symbolic_pauli_exp(&sparse(&[z(0), z(1)]), branch); + }); + let simulation_action = phased_action_from_simulation(&simulation, &[0, 1], &[0, 1]).expect("simulation action"); + + simulation_action + .is_equivalent(&circuit_action) + .expect("simulator-native action must match the circuit action"); +} + +#[test] +fn simulator_native_distinguishes_opposite_signs() { + let positive = choi_simulation(1, |simulation, branch| { + simulation.symbolic_pauli_exp(&sparse(&[z(0)]), branch); + }); + let negative = choi_simulation(1, |simulation, branch| { + simulation.symbolic_pauli_exp(&-sparse(&[z(0)]), branch); + }); + + let positive_action = phased_action_from_simulation(&positive, &[0], &[0]).expect("positive action"); + let negative_action = phased_action_from_simulation(&negative, &[0], &[0]).expect("negative action"); + + positive_action + .is_equivalent_up_to_signs(&negative_action) + .expect("phaseless actions must be identical"); + let reasons = positive_action + .is_equivalent(&negative_action) + .expect_err("opposite signs differ only in relative phase"); + assert_eq!(reasons, vec![ActionsInequivalenceReason::RelativePhase]); +} + +// ================================================================================================ +// "Ejection": remote/measurement-based execution of a Z-diagonal channel. +// +// A Z-diagonal channel on `n` system qubits is applied indirectly: `n` ancillas (prepared in |0⟩) +// receive a transversal CNOT from the system qubits, the Z-diagonal channel acts on the ancillas, +// each ancilla is destructively measured in the X basis, and a "−" outcome triggers a conditional Z +// correction on the corresponding system qubit. The action must equal applying the same Z-diagonal +// channel directly to the system qubits. The X-basis measurements introduce *true* random bits that +// must be marginalized, while the channel's symbolic rotation angles are *virtual* bits that must +// correspond one-to-one — exactly the mixed case the virtual/true distinction is built for. +// ================================================================================================ + +/// A tensor product of `Z` operators, one on each qubit `support[i]` selected by `local_indices`. +/// `local_indices` name positions *within* `support`, letting the same channel description +/// (`angle_supports`) be applied either to the system qubits or to the ancillas. +fn z_product(local_indices: &[usize], support: &[QubitId]) -> SparsePauli { + let positioned: Vec = local_indices.iter().map(|&index| z(support[index])).collect(); + (&positioned[..]).into() +} + +/// Applies the symbolic Z-rotations indexed by `angle_supports` (each a tensor product of `Z`s, with +/// its own symbolic angle) to the qubits named by `support`, in allocation order. +fn apply_symbolic_z_rotations(builder: &mut CircuitBuilder, angle_supports: &[Vec], support: &[QubitId]) { + for local_indices in angle_supports { + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&z_product(local_indices, support), angle); + } +} + +/// The direct circuit: the symbolic Z-rotations applied straight to the `n` system qubits. +fn direct_z_channel(n: usize, angle_supports: &[Vec]) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + apply_symbolic_z_rotations(builder, angle_supports, &system); + }) +} + +/// The ejection circuit: the same Z-rotations are executed *remotely* on `n` fresh ancillas and then +/// teleported back onto the system. Qubits `0..n` are the system, `n..2n` the ancillas. Three phases: +/// 1. entangle — a transversal CNOT copies each system qubit onto its ancilla, +/// 2. rotate remotely — the symbolic Z-rotations act on the ancillas instead of the system, +/// 3. measure and correct — each ancilla is measured in the X basis (a *true* random bit), and a +/// `−` outcome triggers a conditional `Z` correction on the matching system qubit. +/// +/// The net channel on the system must equal [`direct_z_channel`]. +fn z_ejection_channel(n: usize, angle_supports: &[Vec]) -> Circuit { + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + build_circuit(|builder| { + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[system_qubit, ancilla]); + } + apply_symbolic_z_rotations(builder, angle_supports, &ancillas); + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + let minus_outcome = builder.measure(&sparse(&[x(ancilla)])); + builder.conditional_pauli(&sparse(&[z(system_qubit)]), &[minus_outcome], true); + } + }) +} + +/// Asserts the ejection invariant: executing the Z-rotations remotely and teleporting them back is +/// the same channel — including every branch phase — as applying them directly. Checked both ways so +/// the equivalence is symmetric. +fn check_z_ejection(n: usize, angle_supports: &[Vec]) { + let system: Vec = (0..n).collect(); + let direct = direct_z_channel(n, angle_supports); + let ejection = z_ejection_channel(n, angle_supports); + + let direct_action = phased_action_of(&direct, &system, &system).expect("direct channel action"); + let ejection_action = phased_action_of(&ejection, &system, &system).expect("ejection channel action"); + + direct_action.is_equivalent(&ejection_action).unwrap_or_else(|reasons| { + panic!("ejection of {angle_supports:?} on {n} qubits must equal the direct channel: {reasons:?}") + }); + ejection_action + .is_equivalent(&direct_action) + .expect("ejection equivalence must be symmetric"); +} + +#[test] +fn single_qubit_z_rotation_ejection() { + check_z_ejection(1, &[vec![0]]); +} + +#[test] +fn two_qubit_z_rotation_ejections() { + check_z_ejection(2, &[vec![0]]); + check_z_ejection(2, &[vec![1]]); + check_z_ejection(2, &[vec![0, 1]]); + check_z_ejection(2, &[vec![0], vec![1], vec![0, 1]]); +} + +#[test] +fn three_qubit_all_z_products_ejection() { + let all_nontrivial: Vec> = (1u32..8) + .map(|mask| (0..3).filter(|bit| mask & (1 << bit) != 0).collect()) + .collect(); + check_z_ejection(3, &all_nontrivial); +} + +#[test] +fn repeated_angles_ejection() { + check_z_ejection(2, &[vec![0], vec![0], vec![0, 1], vec![0, 1]]); +} + +/// A wrong correction (conditioning on the wrong measurement outcome) must be detected: the branch +/// phase then depends on a true measurement bit that the direct channel cannot reproduce. +#[test] +fn miscorrected_ejection_is_detected() { + let n = 2; + let angle_supports = [vec![0, 1]]; + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + let broken = build_circuit(|builder| { + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[system_qubit, ancilla]); + } + apply_symbolic_z_rotations(builder, &angle_supports, &ancillas); + let mut outcomes = Vec::new(); + for &ancilla in &ancillas { + outcomes.push(builder.measure(&sparse(&[x(ancilla)]))); + } + // Apply only the first correction, dropping the second: leaves a residual outcome dependence. + builder.conditional_pauli(&sparse(&[z(system[0])]), &[outcomes[0]], true); + }); + + let direct = direct_z_channel(n, &angle_supports); + let direct_action = phased_action_of(&direct, &system, &system).expect("direct action"); + let broken_action = phased_action_of(&broken, &system, &system).expect("broken action"); + + direct_action + .is_equivalent(&broken_action) + .expect_err("a missing correction must make the ejection inequivalent"); +} + +// A Z-diagonal *Clifford* (here `S` on each ancilla plus a `CZ`) carries a non-trivial phase but no +// symbolic angles. Ejecting it must equal applying it directly — the **no-angle** case, where the +// phased equivalence reduces exactly to the phaseless `OutcomeCompleteSimulation` behaviour: the only +// random bits are the corrected true X-measurement outcomes, so the relative-phase check is vacuous +// and the residual `|±⟩` ancillas (left uncleaned, exactly as in the phaseless ejection precedent) +// do not affect the comparison. + +fn apply_z_diagonal_clifford(builder: &mut CircuitBuilder, support: &[QubitId]) { + for &qubit in support { + builder.unitary_op(UnitaryOp::SqrtZ, &[qubit]); + } + for window in support.windows(2) { + builder.unitary_op(UnitaryOp::ControlledZ, &[window[0], window[1]]); + } +} + +fn direct_z_clifford_channel(n: usize) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| apply_z_diagonal_clifford(builder, &system)) +} + +fn z_clifford_ejection_channel(n: usize) -> Circuit { + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + build_circuit(|builder| { + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[system_qubit, ancilla]); + } + apply_z_diagonal_clifford(builder, &ancillas); + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + let outcome = builder.measure(&sparse(&[x(ancilla)])); + builder.conditional_pauli(&sparse(&[z(system_qubit)]), &[outcome], true); + } + }) +} + +#[test] +fn z_diagonal_clifford_ejection_without_angles() { + for n in 1..=3 { + let system: Vec = (0..n).collect(); + let direct = direct_z_clifford_channel(n); + let ejection = z_clifford_ejection_channel(n); + let direct_action = phased_action_of(&direct, &system, &system).expect("direct clifford action"); + let ejection_action = phased_action_of(&ejection, &system, &system).expect("ejection clifford action"); + + direct_action.is_equivalent(&ejection_action).unwrap_or_else(|reasons| { + panic!("no-angle Z-diagonal Clifford ejection on {n} qubits must equal direct: {reasons:?}") + }); + ejection_action + .is_equivalent(&direct_action) + .expect("no-angle ejection equivalence must be symmetric"); + } +} + +// ================================================================================================ +// X-basis "ejection": the Hadamard dual of the Z-basis gadget above. +// +// An X-diagonal channel on `n` system qubits is applied indirectly: `n` ancillas (prepared in |+⟩) +// drive a transversal CNOT *into* the system qubits (control = ancilla, target = system), the +// X-diagonal channel acts on the ancillas, each ancilla is destructively measured in the Z basis, +// and a `1` outcome triggers a conditional X correction on the corresponding system qubit. Because +// this whole gadget is the conjugation of the (already verified) Z-basis gadget by a Hadamard on +// every system and ancilla qubit, it must equal applying the same X-diagonal channel directly. As in +// the Z case, the Z-basis measurements introduce *true* random bits that must be marginalized, while +// the symbolic rotation angles are *virtual* bits that correspond one-to-one. +// ================================================================================================ + +/// `X` on each `qubits[i]`-th entry of `support` (a tensor product of `X` operators). +fn x_product(qubits: &[usize], support: &[QubitId]) -> SparsePauli { + let positioned: Vec = qubits.iter().map(|&qubit| x(support[qubit])).collect(); + (&positioned[..]).into() +} + +/// Applies the symbolic X-rotations indexed by `angle_supports` (each a tensor product of `X`s, with +/// its own symbolic angle) to the qubits named by `support`, in allocation order. +fn apply_symbolic_x_rotations(builder: &mut CircuitBuilder, angle_supports: &[Vec], support: &[QubitId]) { + for qubits in angle_supports { + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&x_product(qubits, support), angle); + } +} + +/// The direct circuit: the symbolic X-rotations applied straight to the `n` system qubits. +fn direct_x_channel(n: usize, angle_supports: &[Vec]) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + apply_symbolic_x_rotations(builder, angle_supports, &system); + }) +} + +/// The ejection circuit: the same symbolic X-rotations executed remotely on `n` ancillas. +fn x_ejection_channel(n: usize, angle_supports: &[Vec]) -> Circuit { + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + build_circuit(|builder| { + for &ancilla in &ancillas { + builder.unitary_op(UnitaryOp::Hadamard, &[ancilla]); + } + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[ancilla, system_qubit]); + } + apply_symbolic_x_rotations(builder, angle_supports, &ancillas); + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + let outcome = builder.measure(&sparse(&[z(ancilla)])); + builder.conditional_pauli(&sparse(&[x(system_qubit)]), &[outcome], true); + } + }) +} + +fn check_x_ejection(n: usize, angle_supports: &[Vec]) { + let system: Vec = (0..n).collect(); + let direct = direct_x_channel(n, angle_supports); + let ejection = x_ejection_channel(n, angle_supports); + + let direct_action = phased_action_of(&direct, &system, &system).expect("direct channel action"); + let ejection_action = phased_action_of(&ejection, &system, &system).expect("ejection channel action"); + + direct_action.is_equivalent(&ejection_action).unwrap_or_else(|reasons| { + panic!("X ejection of {angle_supports:?} on {n} qubits must equal the direct channel: {reasons:?}") + }); + ejection_action + .is_equivalent(&direct_action) + .expect("X ejection equivalence must be symmetric"); +} + +#[test] +fn single_qubit_x_rotation_ejection() { + check_x_ejection(1, &[vec![0]]); +} + +#[test] +fn two_qubit_x_rotation_ejections() { + check_x_ejection(2, &[vec![0]]); + check_x_ejection(2, &[vec![1]]); + check_x_ejection(2, &[vec![0, 1]]); + check_x_ejection(2, &[vec![0], vec![1], vec![0, 1]]); +} + +#[test] +fn three_qubit_all_x_products_ejection() { + let all_nontrivial: Vec> = (1u32..8) + .map(|mask| (0..3).filter(|bit| mask & (1 << bit) != 0).collect()) + .collect(); + check_x_ejection(3, &all_nontrivial); +} + +#[test] +fn repeated_x_angles_ejection() { + check_x_ejection(2, &[vec![0], vec![0], vec![0, 1], vec![0, 1]]); +} + +// ================================================================================================ +// Ejection of a *general* Z-diagonal channel: symbolic Z-rotations (virtual bits) mixed with +// non-destructive Z-basis measurements (true observed bits). Ejecting the channel onto ancillas adds +// a third kind of bit, the destructive X-readout outcomes (true auxiliary bits, marginalized). The +// default `is_equivalent` resolves all three automatically: the virtual angle bits correspond one to +// one, the observed measurement bits map identity-by-allocation-order, and the readout bits are +// projected out. This is the first gadget that exercises all three provenance classes at once. +// ================================================================================================ + +/// Applies a Z-diagonal *channel* to `support`: first the symbolic Z-rotations indexed by +/// `angle_supports`, then non-destructive Z-basis measurements of the tensor products indexed by +/// `measure_supports`, all in allocation order. +fn apply_z_diagonal_channel( + builder: &mut CircuitBuilder, + angle_supports: &[Vec], + measure_supports: &[Vec], + support: &[QubitId], +) { + apply_symbolic_z_rotations(builder, angle_supports, support); + for qubits in measure_supports { + let _ = builder.measure(&z_product(qubits, support)); + } +} + +fn direct_z_channel_with_measurements( + n: usize, + angle_supports: &[Vec], + measure_supports: &[Vec], +) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + apply_z_diagonal_channel(builder, angle_supports, measure_supports, &system); + }) +} + +fn z_ejection_channel_with_measurements( + n: usize, + angle_supports: &[Vec], + measure_supports: &[Vec], +) -> Circuit { + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + build_circuit(|builder| { + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[system_qubit, ancilla]); + } + apply_z_diagonal_channel(builder, angle_supports, measure_supports, &ancillas); + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + let outcome = builder.measure(&sparse(&[x(ancilla)])); + builder.conditional_pauli(&sparse(&[z(system_qubit)]), &[outcome], true); + } + }) +} + +fn check_z_ejection_with_measurements(n: usize, angle_supports: &[Vec], measure_supports: &[Vec]) { + let system: Vec = (0..n).collect(); + let direct = direct_z_channel_with_measurements(n, angle_supports, measure_supports); + let ejection = z_ejection_channel_with_measurements(n, angle_supports, measure_supports); + + let direct_action = phased_action_of(&direct, &system, &system).expect("direct channel action"); + let ejection_action = phased_action_of(&ejection, &system, &system).expect("ejection channel action"); + + direct_action.is_equivalent(&ejection_action).unwrap_or_else(|reasons| { + panic!("Z ejection of angles {angle_supports:?} and measurements {measure_supports:?} on {n} qubits must equal the direct channel: {reasons:?}") + }); + ejection_action + .is_equivalent(&direct_action) + .expect("Z channel ejection equivalence must be symmetric"); +} + +#[test] +fn single_qubit_z_channel_ejection() { + check_z_ejection_with_measurements(1, &[vec![0]], &[vec![0]]); + check_z_ejection_with_measurements(1, &[], &[vec![0]]); +} + +#[test] +fn two_qubit_z_channel_ejections() { + check_z_ejection_with_measurements(2, &[vec![0]], &[vec![1]]); + check_z_ejection_with_measurements(2, &[vec![0], vec![1]], &[vec![0, 1]]); + check_z_ejection_with_measurements(2, &[vec![0, 1]], &[vec![0], vec![1]]); + check_z_ejection_with_measurements(2, &[], &[vec![0], vec![1], vec![0, 1]]); +} + +#[test] +fn three_qubit_z_channel_ejection() { + check_z_ejection_with_measurements(3, &[vec![0, 1, 2]], &[vec![0], vec![1, 2]]); +} + +// ================================================================================================ +// X-basis dual of the general-channel ejection above: symbolic X-rotations mixed with +// non-destructive X-basis measurements, ejected through `|+⟩` ancillas with reversed CNOTs, +// destructive Z-readout, and conditional X corrections. +// ================================================================================================ + +/// Applies an X-diagonal *channel* to `support`: the symbolic X-rotations indexed by +/// `angle_supports`, then non-destructive X-basis measurements indexed by `measure_supports`. +fn apply_x_diagonal_channel( + builder: &mut CircuitBuilder, + angle_supports: &[Vec], + measure_supports: &[Vec], + support: &[QubitId], +) { + apply_symbolic_x_rotations(builder, angle_supports, support); + for qubits in measure_supports { + let _ = builder.measure(&x_product(qubits, support)); + } +} + +fn direct_x_channel_with_measurements( + n: usize, + angle_supports: &[Vec], + measure_supports: &[Vec], +) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + apply_x_diagonal_channel(builder, angle_supports, measure_supports, &system); + }) +} + +fn x_ejection_channel_with_measurements( + n: usize, + angle_supports: &[Vec], + measure_supports: &[Vec], +) -> Circuit { + let system: Vec = (0..n).collect(); + let ancillas: Vec = (n..2 * n).collect(); + build_circuit(|builder| { + for &ancilla in &ancillas { + builder.unitary_op(UnitaryOp::Hadamard, &[ancilla]); + } + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + builder.unitary_op(UnitaryOp::ControlledX, &[ancilla, system_qubit]); + } + apply_x_diagonal_channel(builder, angle_supports, measure_supports, &ancillas); + for (&system_qubit, &ancilla) in system.iter().zip(ancillas.iter()) { + let outcome = builder.measure(&sparse(&[z(ancilla)])); + builder.conditional_pauli(&sparse(&[x(system_qubit)]), &[outcome], true); + } + }) +} + +fn check_x_ejection_with_measurements(n: usize, angle_supports: &[Vec], measure_supports: &[Vec]) { + let system: Vec = (0..n).collect(); + let direct = direct_x_channel_with_measurements(n, angle_supports, measure_supports); + let ejection = x_ejection_channel_with_measurements(n, angle_supports, measure_supports); + + let direct_action = phased_action_of(&direct, &system, &system).expect("direct channel action"); + let ejection_action = phased_action_of(&ejection, &system, &system).expect("ejection channel action"); + + direct_action.is_equivalent(&ejection_action).unwrap_or_else(|reasons| { + panic!("X ejection of angles {angle_supports:?} and measurements {measure_supports:?} on {n} qubits must equal the direct channel: {reasons:?}") + }); + ejection_action + .is_equivalent(&direct_action) + .expect("X channel ejection equivalence must be symmetric"); +} + +#[test] +fn single_qubit_x_channel_ejection() { + check_x_ejection_with_measurements(1, &[vec![0]], &[vec![0]]); + check_x_ejection_with_measurements(1, &[], &[vec![0]]); +} + +#[test] +fn two_qubit_x_channel_ejections() { + check_x_ejection_with_measurements(2, &[vec![0]], &[vec![1]]); + check_x_ejection_with_measurements(2, &[vec![0], vec![1]], &[vec![0, 1]]); + check_x_ejection_with_measurements(2, &[vec![0, 1]], &[vec![0], vec![1]]); + check_x_ejection_with_measurements(2, &[], &[vec![0], vec![1], vec![0, 1]]); +} + +#[test] +fn three_qubit_x_channel_ejection() { + check_x_ejection_with_measurements(3, &[vec![0, 1, 2]], &[vec![0], vec![1, 2]]); +} + +// ================================================================================================ +// Section 4.1 of arXiv:2603.24717: verifying parameterized state-preparation circuits. +// +// To decide whether two parameterized circuits prepare the same state for *every* rotation angle, +// C₁ exp(iα Z) C₂|0…0> == D₁ exp(iα Z) D₂|0…0> (for all α), +// it suffices to check a single EXACT stabilizer-state equality with the angle replaced by a binary +// symbolic exponent, +// C₁ Z^a C₂|0…0> == D₁ Z^a D₂|0…0>, +// because exactness — equality including the relative phase between the a = 0 and a = 1 branches — +// pins down the rotation phase for every α. This needs no dedicated verification entry point: the +// check is exactly `phased_action_of` + `PhasedCircuitAction::is_equivalent`, the phased analog of +// how `OutcomeCompleteSimulation` performs phaseless equality checking. `Z^a` is realized by an +// `allocate_symbolic_angle` angle applied with `symbolic_pauli_exp`. +// ================================================================================================ + +/// Records a state-preparation gadget `C₁ (∏ₖ Z_k^{a_k}) C₂ |0…0>` as a phased action with no input +/// qubits — the `inputs = []` (state-preparation) case of `phased_action_of`. +fn prepared_state_action(qubit_count: usize, build: impl FnOnce(&mut CircuitBuilder)) -> PhasedCircuitAction { + let outputs: Vec = (0..qubit_count).collect(); + let circuit = build_circuit(build); + phased_action_of(&circuit, &[], &outputs).expect("state preparation action") +} + +/// Prepares `|+…+>` by Hadamarding every qubit in `0..qubit_count`. +fn prepare_plus(builder: &mut CircuitBuilder, qubit_count: usize) { + for qubit in 0..qubit_count { + builder.unitary_op(UnitaryOp::Hadamard, &[qubit]); + } +} + +/// Two different Clifford factorizations `C₁ Z^a C₂` of the same parameterized state must verify as +/// equivalent: `exp(iα Z₀Z₁)|++>` realized directly, versus the CNOT-conjugated single-qubit rotation +/// `CNOT₀₁ exp(iα Z₁) CNOT₀₁ |++>` (using `CNOT₀₁ Z₁ CNOT₀₁ = Z₀Z₁` and `CNOT₀₁|++> = |++>`). +#[test] +fn verifies_equal_state_preparation_factorizations() { + let direct = prepared_state_action(2, |builder| { + prepare_plus(builder, 2); + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0), z(1)]), angle); + }); + let conjugated = prepared_state_action(2, |builder| { + prepare_plus(builder, 2); + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(1)]), angle); + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + }); + + direct + .is_equivalent(&conjugated) + .expect("§4.1: the two factorizations prepare the same parameterized state"); + conjugated + .is_equivalent(&direct) + .expect("verification must be symmetric"); +} + +/// The check is phase-sensitive: `exp(+iα Z₀)|+>` and `exp(-iα Z₀)|+>` have identical stabilizer data +/// but opposite branch phase, so they must be distinguished — by exactly one `RelativePhase` reason. +#[test] +fn detects_phase_only_state_preparation_difference() { + let positive = prepared_state_action(1, |builder| { + prepare_plus(builder, 1); + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0)]), angle); + }); + let negative = prepared_state_action(1, |builder| { + prepare_plus(builder, 1); + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&-sparse(&[z(0)]), angle); + }); + + positive + .is_equivalent_up_to_signs(&negative) + .expect("the phaseless data is identical"); + let reasons = positive + .is_equivalent(&negative) + .expect_err("exp(+iαZ) and exp(-iαZ) prepare states differing only in branch phase"); + assert_eq!(reasons, vec![ActionsInequivalenceReason::RelativePhase]); +} + +/// The §4.1 reduction generalizes to several independent symbolic angles. `exp(iα Z₀Z₁) exp(iβ Z₀)|++>` +/// verifies equal to its CNOT-conjugated factorization (angles allocated in the same order, so the +/// virtual-angle bits correspond one to one), while negating the second rotation's Pauli yields a +/// pure branch-phase difference that is detected. +#[test] +fn verifies_multi_angle_state_preparation() { + let direct = |negate_second: bool| { + prepared_state_action(2, move |builder| { + prepare_plus(builder, 2); + let first = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0), z(1)]), first); + let second = builder.allocate_symbolic_angle(); + let pauli = if negate_second { + -sparse(&[z(0)]) + } else { + sparse(&[z(0)]) + }; + builder.symbolic_pauli_exp(&pauli, second); + }) + }; + let conjugated = prepared_state_action(2, |builder| { + prepare_plus(builder, 2); + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + let first = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(1)]), first); + builder.unitary_op(UnitaryOp::ControlledX, &[0, 1]); + let second = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&sparse(&[z(0)]), second); + }); + + direct(false) + .is_equivalent(&conjugated) + .expect("§4.1: the multi-angle factorizations prepare the same parameterized state"); + let reasons = direct(true) + .is_equivalent(&conjugated) + .expect_err("negating one rotation must produce a detectable branch-phase difference"); + assert_eq!(reasons, vec![ActionsInequivalenceReason::RelativePhase]); +} + +// ================================================================================================ +// Negative and randomized tests: a phased equivalence is exact, so any sign flip on a symbolic Pauli +// exponent or any permutation of the symbolic angles on the right-hand side must be detected. The +// phaseless action is unchanged (the symplectic content is identical), so the failure is a pure +// `RelativePhase` (sign flips) or a support/phase mismatch (permutations). Random instances exercise +// the same invariants over many angle supports, qubit counts, and qubit assignments. +// ================================================================================================ + +/// Direct Z-channel, but the `k`-th symbolic Z-rotation is `exp(±iα Z…)` according to `signs[k]`. +fn signed_z_channel(n: usize, angle_supports: &[Vec], signs: &[bool]) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + for (qubits, &negate) in angle_supports.iter().zip(signs.iter()) { + let angle = builder.allocate_symbolic_angle(); + let pauli = if negate { + -z_product(qubits, &system) + } else { + z_product(qubits, &system) + }; + builder.symbolic_pauli_exp(&pauli, angle); + } + }) +} + +/// Direct Z-channel whose `k`-th allocated angle drives the rotation on `angle_supports[perm[k]]`. +fn permuted_z_channel(n: usize, angle_supports: &[Vec], perm: &[usize]) -> Circuit { + let system: Vec = (0..n).collect(); + build_circuit(|builder| { + for &target in perm { + let angle = builder.allocate_symbolic_angle(); + builder.symbolic_pauli_exp(&z_product(&angle_supports[target], &system), angle); + } + }) +} + +/// All non-trivial Z products on `n` qubits, in ascending-mask order: distinct, independent supports. +fn distinct_z_supports(n: usize) -> Vec> { + (1u32..(1 << n)) + .map(|mask| (0..n).filter(|bit| mask & (1 << bit) != 0).collect()) + .collect() +} + +#[test] +fn flipping_any_sign_yields_relative_phase() { + let n = 2; + let supports = distinct_z_supports(n); + let system: Vec = (0..n).collect(); + let baseline = signed_z_channel(n, &supports, &vec![false; supports.len()]); + let baseline_action = phased_action_of(&baseline, &system, &system).expect("baseline action"); + for mask in 1u32..(1 << supports.len()) { + let signs: Vec = (0..supports.len()).map(|k| mask & (1 << k) != 0).collect(); + let flipped = signed_z_channel(n, &supports, &signs); + let flipped_action = phased_action_of(&flipped, &system, &system).expect("flipped action"); + baseline_action + .is_equivalent_up_to_signs(&flipped_action) + .expect("sign flips leave the phaseless action unchanged"); + let reasons = baseline_action + .is_equivalent(&flipped_action) + .expect_err("sign mask must be detected"); + assert_eq!( + reasons, + vec![ActionsInequivalenceReason::RelativePhase], + "mask {mask:b}" + ); + } +} + +#[test] +fn permuting_distinct_angles_is_detected() { + let n = 2; + let supports = distinct_z_supports(n); + let system: Vec = (0..n).collect(); + let identity: Vec = (0..supports.len()).collect(); + let base = permuted_z_channel(n, &supports, &identity); + let base_action = phased_action_of(&base, &system, &system).expect("identity action"); + base_action + .is_equivalent(&base_action) + .expect("identity permutation is self-equivalent"); + let swapped = [1usize, 0, 2]; + let swapped_action = phased_action_of(&permuted_z_channel(n, &supports, &swapped), &system, &system).expect("swap"); + base_action + .is_equivalent(&swapped_action) + .expect_err("permuting distinct angles must be inequivalent"); +} + +prop_compose! { + fn arbitrary_signed_z_channel(qubit_range: Range) + (n in qubit_range)(signs in proptest::collection::vec(any::(), 1..(1usize << n)), n in Just(n)) + -> (usize, Vec) + { (n, signs) } +} + +proptest! { + #[test] + fn random_sign_flip_is_pure_relative_phase((n, signs) in arbitrary_signed_z_channel(2..4usize)) { + let supports = distinct_z_supports(n); + let mut signs: Vec = signs.into_iter().take(supports.len()).collect(); + signs.resize(supports.len(), false); + let system: Vec = (0..n).collect(); + let baseline = signed_z_channel(n, &supports, &vec![false; supports.len()]); + let flipped = signed_z_channel(n, &supports, &signs); + let a = phased_action_of(&baseline, &system, &system).expect("baseline"); + let b = phased_action_of(&flipped, &system, &system).expect("flipped"); + a.is_equivalent_up_to_signs(&b).expect("phaseless equal"); + match a.is_equivalent(&b) { + Ok(()) => prop_assert!(signs.iter().all(|s| !s), "only the all-positive mask is equivalent"), + Err(reasons) => prop_assert_eq!(reasons, vec![ActionsInequivalenceReason::RelativePhase]), + } + } +} + +prop_compose! { + fn arbitrary_permutation(qubit_range: Range)(n in qubit_range, seed in any::()) -> (usize, Vec) { + use rand::seq::SliceRandom; + let mut perm: Vec = (0..((1usize << n) - 1)).collect(); + perm.shuffle(&mut rand::rngs::StdRng::seed_from_u64(seed)); + (n, perm) + } +} + +proptest! { + #[test] + fn random_angle_permutation_matches_iff_identity((n, perm) in arbitrary_permutation(2..4usize)) { + let supports = distinct_z_supports(n); + let system: Vec = (0..n).collect(); + let identity: Vec = (0..supports.len()).collect(); + let base = phased_action_of(&permuted_z_channel(n, &supports, &identity), &system, &system).expect("base"); + let permuted = phased_action_of(&permuted_z_channel(n, &supports, &perm), &system, &system).expect("perm"); + if perm == identity { + prop_assert!(base.is_equivalent(&permuted).is_ok()); + } else { + prop_assert!(base.is_equivalent(&permuted).is_err(), "permutation {perm:?} maps distinct angles to wrong Paulis"); + } + } + + #[test] + fn random_multi_angle_channel_self_equivalent((n, signs) in arbitrary_signed_z_channel(2..4usize)) { + let supports = distinct_z_supports(n); + let mut signs: Vec = signs.into_iter().take(supports.len()).collect(); + signs.resize(supports.len(), false); + let system: Vec = (0..n).collect(); + let a = phased_action_of(&signed_z_channel(n, &supports, &signs), &system, &system).expect("a"); + let b = phased_action_of(&signed_z_channel(n, &supports, &signs), &system, &system).expect("b"); + prop_assert!(a.is_equivalent(&b).is_ok(), "a channel must be exactly equivalent to itself"); + } +} diff --git a/pauliverse/tests/phased_outcome_complete_dense.rs b/pauliverse/tests/phased_outcome_complete_dense.rs new file mode 100644 index 00000000..3c2bf0ca --- /dev/null +++ b/pauliverse/tests/phased_outcome_complete_dense.rs @@ -0,0 +1,480 @@ +//! Dense statevector validation of `PhasedOutcomeCompleteSimulation`. +//! +//! For random circuits (Clifford gates, Pauli exponentials, Paulis, controlled-Paulis, conditional +//! Paulis and Pauli measurements) this test enumerates every random-bit assignment `r`, materialises +//! the simulator's claimed output state `i^⟨p,r⟩ (-1)^⟨Br+s,r⟩ R|Ar⟩`, and compares it — *exactly, +//! including the global phase* — against a brute-force dense statevector simulation whose +//! measurement outcomes are forced to the branch selected by `r`. + +#![allow( + clippy::cast_precision_loss, + clippy::cast_possible_truncation, + clippy::cast_sign_loss, + clippy::too_many_lines +)] + +use binar::vec::AlignedBitVec; +use binar::{Bitwise, BitwiseMut}; +use paulimer::clifford::Clifford; +use paulimer::pauli::commutes_with; +use paulimer::{DensePauli, SparsePauli, UnitaryOp}; +use pauliverse::{PhasedOutcomeCompleteSimulation, Simulation}; +use rand::RngExt; + +use dense_oracle::{C, Dense, close, gate_matrix, normalize, pauli_arrays, statevector, zeta8}; + +#[derive(Clone)] +enum Op { + Gate(UnitaryOp, Vec), + Pauli(String), + PauliExp(String), + ControlledPauli(String, String), + ConditionalPauli(String, Vec, bool), + Measure(String), +} + +fn random_hermitian_pauli(rng: &mut impl RngExt, qubit_count: usize) -> String { + loop { + let mut letters = String::new(); + let mut any = false; + for _ in 0..qubit_count { + match rng.random_range(0..4) { + 0 => letters.push('I'), + 1 => { + letters.push('X'); + any = true; + } + 2 => { + letters.push('Z'); + any = true; + } + _ => { + letters.push('Y'); + any = true; + } + } + } + if !any { + continue; + } + let sign = if rng.random_range(0..2) == 0 { "-" } else { "" }; + return format!("{sign}{letters}"); + } +} + +fn two_distinct(rng: &mut impl RngExt, qubit_count: usize) -> (usize, usize) { + let first = rng.random_range(0..qubit_count); + let mut second = rng.random_range(0..qubit_count); + while second == first { + second = rng.random_range(0..qubit_count); + } + (first, second) +} + +fn random_circuit(rng: &mut impl RngExt, qubit_count: usize) -> Vec { + let single_qubit_gates = [ + UnitaryOp::Hadamard, + UnitaryOp::X, + UnitaryOp::Y, + UnitaryOp::Z, + UnitaryOp::SqrtZ, + UnitaryOp::SqrtZInv, + UnitaryOp::SqrtX, + UnitaryOp::SqrtXInv, + UnitaryOp::SqrtY, + UnitaryOp::SqrtYInv, + ]; + let two_qubit_gates = [UnitaryOp::ControlledX, UnitaryOp::ControlledZ, UnitaryOp::Swap]; + let mut ops = Vec::new(); + let mut measurement_count = 0usize; + let op_count = rng.random_range(6..14); + for _ in 0..op_count { + match rng.random_range(0..7) { + 0 => { + let qubit = rng.random_range(0..qubit_count); + ops.push(Op::Gate( + single_qubit_gates[rng.random_range(0..single_qubit_gates.len())], + vec![qubit], + )); + } + 1 => { + let (first_qubit, second_qubit) = two_distinct(rng, qubit_count); + ops.push(Op::Gate( + two_qubit_gates[rng.random_range(0..two_qubit_gates.len())], + vec![first_qubit, second_qubit], + )); + } + 2 => ops.push(Op::Pauli(random_hermitian_pauli(rng, qubit_count))), + 3 => ops.push(Op::PauliExp(random_hermitian_pauli(rng, qubit_count))), + 4 => { + let first_pauli_string = random_hermitian_pauli(rng, qubit_count); + let mut second_pauli_string = random_hermitian_pauli(rng, qubit_count); + let mut guard = 0; + loop { + let first_sparse_pauli: SparsePauli = first_pauli_string.parse().unwrap(); + let second_sparse_pauli: SparsePauli = second_pauli_string.parse().unwrap(); + if commutes_with(&first_sparse_pauli, &second_sparse_pauli) { + break; + } + second_pauli_string = random_hermitian_pauli(rng, qubit_count); + guard += 1; + if guard > 32 { + break; + } + } + let first_sparse_pauli: SparsePauli = first_pauli_string.parse().unwrap(); + let second_sparse_pauli: SparsePauli = second_pauli_string.parse().unwrap(); + if commutes_with(&first_sparse_pauli, &second_sparse_pauli) { + ops.push(Op::ControlledPauli(first_pauli_string, second_pauli_string)); + } + } + 5 => { + if measurement_count > 0 && rng.random_range(0..2) == 0 { + let mut outcomes = Vec::new(); + for outcome_index in 0..measurement_count { + if rng.random_range(0..2) == 0 { + outcomes.push(outcome_index); + } + } + if !outcomes.is_empty() { + ops.push(Op::ConditionalPauli( + random_hermitian_pauli(rng, qubit_count), + outcomes, + rng.random_range(0..2) == 1, + )); + } + } + } + _ => { + if measurement_count < 5 { + ops.push(Op::Measure(random_hermitian_pauli(rng, qubit_count))); + measurement_count += 1; + } + } + } + } + ops +} + +fn run_simulation(ops: &[Op], qubit_count: usize) -> PhasedOutcomeCompleteSimulation { + let mut sim = PhasedOutcomeCompleteSimulation::new(qubit_count); + for op in ops { + match op { + Op::Gate(gate, support) => sim.unitary_op(*gate, support), + Op::Pauli(pauli) => sim.pauli(&pauli.parse().unwrap()), + Op::PauliExp(pauli) => sim.pauli_exp(&pauli.parse().unwrap()), + Op::ControlledPauli(first_pauli, second_pauli) => { + sim.controlled_pauli(&first_pauli.parse().unwrap(), &second_pauli.parse().unwrap()); + } + Op::ConditionalPauli(pauli, outcomes, parity) => { + sim.conditional_pauli(&pauli.parse().unwrap(), outcomes, *parity); + } + Op::Measure(pauli) => { + sim.measure(&pauli.parse().unwrap()); + } + } + } + sim +} + +fn dense_reference(ops: &[Op], outcome_bits: &[bool], qubit_count: usize) -> Vec { + let mut dense = Dense::zero(qubit_count); + let mut measurement_index = 0usize; + for op in ops { + match op { + Op::Gate(gate, support) => match gate { + UnitaryOp::ControlledX => dense.apply_cx(support[0], support[1]), + UnitaryOp::ControlledZ => dense.apply_cz(support[0], support[1]), + UnitaryOp::Swap => dense.apply_swap(support[0], support[1]), + other => dense.apply1(support[0], gate_matrix(*other)), + }, + Op::Pauli(pauli) => { + let (x_bits, z_bits, phase) = pauli_arrays(&pauli.parse::().unwrap(), qubit_count); + dense.apply_pauli(&x_bits, &z_bits, phase); + } + Op::PauliExp(pauli) => { + let (x_bits, z_bits, phase) = pauli_arrays(&pauli.parse::().unwrap(), qubit_count); + dense.apply_pauli_exp(&x_bits, &z_bits, phase); + } + Op::ControlledPauli(first_pauli, second_pauli) => { + let first_arrays = pauli_arrays(&first_pauli.parse::().unwrap(), qubit_count); + let second_arrays = pauli_arrays(&second_pauli.parse::().unwrap(), qubit_count); + dense.apply_controlled_pauli(&first_arrays, &second_arrays); + } + Op::ConditionalPauli(pauli, outcomes, parity) => { + let condition = outcomes + .iter() + .fold(false, |acc, &outcome_index| acc ^ outcome_bits[outcome_index]); + if condition == *parity { + let (x_bits, z_bits, phase) = pauli_arrays(&pauli.parse::().unwrap(), qubit_count); + dense.apply_pauli(&x_bits, &z_bits, phase); + } + } + Op::Measure(pauli) => { + let (x_bits, z_bits, phase) = pauli_arrays(&pauli.parse::().unwrap(), qubit_count); + dense.project(&x_bits, &z_bits, phase, outcome_bits[measurement_index]); + measurement_index += 1; + } + } + } + dense.amp +} + +fn claimed_state(sim: &PhasedOutcomeCompleteSimulation, random_bits: &[bool], qubit_count: usize) -> Vec { + let encoder = sim.phased_state_encoder(); + let base = statevector(&encoder); + + let sign_matrix = sim.aligned_sign_matrix(); + let random_outcome_count = sim.random_outcome_count(); + let mut register = AlignedBitVec::zeros(qubit_count); + for qubit in 0..qubit_count { + let mut bit = false; + for (column, random_bit) in random_bits.iter().enumerate().take(random_outcome_count) { + if *random_bit && sign_matrix[(qubit, column)] { + bit = !bit; + } + } + register.assign_index(qubit, bit); + } + + let image = encoder.clifford().image_x_bits(®ister); + let (x_bits, z_bits, phase) = pauli_arrays(&image, qubit_count); + let mut dense = Dense { qubit_count, amp: base }; + dense.apply_pauli(&x_bits, &z_bits, phase); + + let exponent = i64::from(sim.output_phase_exponent(random_bits)); + for amplitude in &mut dense.amp { + *amplitude *= zeta8(exponent); + } + let mut amp = dense.amp; + normalize(&mut amp); + amp +} + +fn outcome_vector(sim: &PhasedOutcomeCompleteSimulation, random_bits: &[bool]) -> Vec { + let outcome_matrix = sim.aligned_outcome_matrix(); + let shift = sim.aligned_outcome_shift(); + let random_outcome_count = sim.random_outcome_count(); + (0..sim.outcome_count()) + .map(|row| { + let mut bit = shift.index(row); + for (column, random_bit) in random_bits.iter().enumerate().take(random_outcome_count) { + if *random_bit && outcome_matrix.row(row).index(column) { + bit = !bit; + } + } + bit + }) + .collect() +} + +fn describe(ops: &[Op]) -> String { + ops.iter() + .map(|op| match op { + Op::Gate(gate, support) => format!("Gate({gate:?},{support:?})"), + Op::Pauli(pauli) => format!("Pauli({pauli})"), + Op::PauliExp(pauli) => format!("PauliExp({pauli})"), + Op::ControlledPauli(first_pauli, second_pauli) => format!("CPauli({first_pauli},{second_pauli})"), + Op::ConditionalPauli(pauli, outcomes, parity) => format!("CondPauli({pauli},{outcomes:?},{parity})"), + Op::Measure(pauli) => format!("Measure({pauli})"), + }) + .collect::>() + .join(" | ") +} + +fn verify(ops: &[Op], qubit_count: usize) { + let sim = run_simulation(ops, qubit_count); + let random_outcome_count = sim.random_outcome_count(); + assert!(random_outcome_count <= 12, "too many random bits to enumerate"); + for assignment in 0..(1usize << random_outcome_count) { + let random_bits: Vec = (0..random_outcome_count) + .map(|bit| (assignment >> bit) & 1 == 1) + .collect(); + let outcomes = outcome_vector(&sim, &random_bits); + let reference = dense_reference(ops, &outcomes, qubit_count); + let claimed = claimed_state(&sim, &random_bits, qubit_count); + assert!( + close(&claimed, &reference), + "mismatch: ops=[{}] random_bits={random_bits:?}", + describe(ops) + ); + } +} + +#[test] +fn single_pauli_measurement() { + verify(&[Op::Measure("X".into())], 1); + verify(&[Op::Measure("Y".into())], 1); + verify(&[Op::Measure("-X".into())], 1); + verify(&[Op::Measure("-Y".into())], 1); +} + +#[test] +fn two_pauli_measurements() { + verify(&[Op::Measure("X".into()), Op::Measure("Y".into())], 1); + verify(&[Op::Measure("Y".into()), Op::Measure("X".into())], 1); + verify(&[Op::Measure("X".into()), Op::Measure("Z".into())], 1); + verify(&[Op::Measure("-X".into()), Op::Measure("-Y".into())], 1); +} + +#[test] +fn measurement_then_conditional() { + verify( + &[ + Op::Measure("X".into()), + Op::ConditionalPauli("Z".into(), vec![0], false), + ], + 1, + ); + verify( + &[Op::Measure("Y".into()), Op::ConditionalPauli("X".into(), vec![0], true)], + 1, + ); +} + +#[test] +fn controlled_pauli_no_randomness() { + verify( + &[ + Op::Gate(UnitaryOp::Hadamard, vec![0]), + Op::ControlledPauli("ZI".into(), "IX".into()), + ], + 2, + ); + verify( + &[ + Op::Gate(UnitaryOp::Hadamard, vec![0]), + Op::ControlledPauli("ZZ".into(), "XX".into()), + ], + 2, + ); + verify( + &[ + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::ControlledPauli("YI".into(), "IY".into()), + ], + 2, + ); +} + +#[test] +fn entangling_then_two_measurements() { + verify( + &[ + Op::Gate(UnitaryOp::Hadamard, vec![0]), + Op::Gate(UnitaryOp::ControlledX, vec![0, 1]), + Op::Measure("XX".into()), + Op::Measure("ZI".into()), + ], + 2, + ); + verify( + &[ + Op::Measure("X".into()), + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::Measure("X".into()), + ], + 1, + ); +} + +#[test] +fn shrink_a() { + verify( + &[ + Op::Gate(UnitaryOp::Swap, vec![2, 0]), + Op::Measure("-ZIX".into()), + Op::Gate(UnitaryOp::ControlledX, vec![2, 1]), + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::ConditionalPauli("-ZXI".into(), vec![0], false), + Op::Gate(UnitaryOp::ControlledZ, vec![2, 1]), + Op::Measure("XIY".into()), + ], + 3, + ); +} + +#[test] +fn shrink_b_no_conditional() { + verify( + &[ + Op::Gate(UnitaryOp::Swap, vec![2, 0]), + Op::Measure("-ZIX".into()), + Op::Gate(UnitaryOp::ControlledX, vec![2, 1]), + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::Gate(UnitaryOp::ControlledZ, vec![2, 1]), + Op::Measure("XIY".into()), + ], + 3, + ); +} + +#[test] +fn shrink_c_conditional_between() { + verify( + &[ + Op::Measure("-ZIX".into()), + Op::ConditionalPauli("-ZXI".into(), vec![0], false), + Op::Measure("XIY".into()), + ], + 3, + ); +} + +#[test] +fn shrink_d_two_meas_gate() { + verify( + &[ + Op::Measure("ZIX".into()), + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::Measure("XIY".into()), + ], + 3, + ); +} + +#[test] +fn captured_regression_one() { + verify( + &[ + Op::Gate(UnitaryOp::Swap, vec![2, 0]), + Op::Measure("-ZIX".into()), + Op::Gate(UnitaryOp::ControlledX, vec![2, 1]), + Op::Gate(UnitaryOp::SqrtX, vec![0]), + Op::ConditionalPauli("-ZXI".into(), vec![0], false), + Op::Gate(UnitaryOp::ControlledZ, vec![2, 1]), + Op::Measure("XIY".into()), + Op::Gate(UnitaryOp::SqrtY, vec![1]), + Op::Gate(UnitaryOp::ControlledX, vec![1, 2]), + Op::Gate(UnitaryOp::Y, vec![2]), + ], + 3, + ); +} + +#[test] +fn phased_outcome_complete_tracks_dense_statevector() { + let mut rng = rand::rng(); + for _trial in 0..600 { + let qubit_count = 3usize; + let ops = random_circuit(&mut rng, qubit_count); + let sim = run_simulation(&ops, qubit_count); + let random_outcome_count = sim.random_outcome_count(); + if random_outcome_count > 8 { + continue; + } + for assignment in 0..(1usize << random_outcome_count) { + let random_bits: Vec = (0..random_outcome_count) + .map(|bit| (assignment >> bit) & 1 == 1) + .collect(); + let outcomes = outcome_vector(&sim, &random_bits); + let reference = dense_reference(&ops, &outcomes, qubit_count); + let claimed = claimed_state(&sim, &random_bits, qubit_count); + assert!( + close(&claimed, &reference), + "mismatch: ops=[{}] random_bits={random_bits:?}", + describe(&ops) + ); + } + } +} diff --git a/test-utils/dense-oracle/Cargo.toml b/test-utils/dense-oracle/Cargo.toml new file mode 100644 index 00000000..b1248cf3 --- /dev/null +++ b/test-utils/dense-oracle/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "dense-oracle" +version = "0.0.0" +edition = "2024" +publish = false +license = "MIT" +description = "Dense (full state vector) simulation oracle shared by paulimer and pauliverse tests" + +[dependencies] +binar = { path = "../../binar", version = "0.1.2" } +paulimer = { path = "../../paulimer", version = "0.2.2" } +num-complex = "0.4" diff --git a/test-utils/dense-oracle/src/lib.rs b/test-utils/dense-oracle/src/lib.rs new file mode 100644 index 00000000..ddd688b9 --- /dev/null +++ b/test-utils/dense-oracle/src/lib.rs @@ -0,0 +1,259 @@ +//! Dense (full state vector) simulation oracle shared by the `paulimer` and +//! `pauliverse` integration tests. +//! +//! This crate exists only to remove the duplicated brute-force statevector +//! simulator that both test suites used to carry inline. It is never published +//! (`publish = false`) and is depended on only as a `dev-dependency`. + +#![allow( + clippy::cast_precision_loss, + clippy::cast_possible_truncation, + clippy::cast_sign_loss, + clippy::too_many_lines, + clippy::needless_range_loop, + clippy::should_implement_trait, + clippy::must_use_candidate, + clippy::return_self_not_must_use +)] + +use binar::{Bitwise, BitwiseMut}; +use num_complex::Complex; +use paulimer::clifford::{Clifford, PhasedCliffordUnitary}; +use paulimer::pauli::Pauli; +use paulimer::{DensePauli, UnitaryOp}; + +/// Complex amplitude type used throughout the dense oracle. +pub type C = Complex; + +/// `exp(i * pi/4 * k)`, i.e. the `k`-th power of the primitive 8th root of unity. +pub fn zeta8(k: i64) -> C { + Complex::cis(std::f64::consts::FRAC_PI_4 * k.rem_euclid(8) as f64) +} + +/// `1 / sqrt(2)`. +pub const ROOT_HALF: f64 = std::f64::consts::FRAC_1_SQRT_2; + +/// Dense amplitude-vector state used as a correctness oracle. +pub struct Dense { + pub qubit_count: usize, + pub amp: Vec, +} + +impl Dense { + pub fn zero(qubit_count: usize) -> Dense { + let mut amp = vec![C::ZERO; 1 << qubit_count]; + amp[0] = C::new(1.0, 0.0); + Dense { qubit_count, amp } + } + pub fn apply1(&mut self, qubit: usize, matrix: [[C; 2]; 2]) { + let bit = 1usize << (self.qubit_count - 1 - qubit); + for base in 0..(1 << self.qubit_count) { + if base & bit == 0 { + let amplitude_0 = self.amp[base]; + let amplitude_1 = self.amp[base | bit]; + self.amp[base] = matrix[0][0] * amplitude_0 + matrix[0][1] * amplitude_1; + self.amp[base | bit] = matrix[1][0] * amplitude_0 + matrix[1][1] * amplitude_1; + } + } + } + pub fn apply_cx(&mut self, control: usize, target: usize) { + let control_bit = 1usize << (self.qubit_count - 1 - control); + let target_bit = 1usize << (self.qubit_count - 1 - target); + let mut out = self.amp.clone(); + for base in 0..(1 << self.qubit_count) { + let src = if base & control_bit != 0 { + base ^ target_bit + } else { + base + }; + out[base] = self.amp[src]; + } + self.amp = out; + } + pub fn apply_cz(&mut self, first_qubit: usize, second_qubit: usize) { + let first_bit = 1usize << (self.qubit_count - 1 - first_qubit); + let second_bit = 1usize << (self.qubit_count - 1 - second_qubit); + for base in 0..(1 << self.qubit_count) { + if base & first_bit != 0 && base & second_bit != 0 { + self.amp[base] = -self.amp[base]; + } + } + } + pub fn apply_swap(&mut self, first_qubit: usize, second_qubit: usize) { + let first_bit = 1usize << (self.qubit_count - 1 - first_qubit); + let second_bit = 1usize << (self.qubit_count - 1 - second_qubit); + let mut out = self.amp.clone(); + for base in 0..(1 << self.qubit_count) { + let bit_first = usize::from(base & first_bit != 0); + let bit_second = usize::from(base & second_bit != 0); + let mut src = base & !first_bit & !second_bit; + if bit_second != 0 { + src |= first_bit; + } + if bit_first != 0 { + src |= second_bit; + } + out[base] = self.amp[src]; + } + self.amp = out; + } + pub fn pauli_applied(&self, x_bits: &[bool], z_bits: &[bool], phase: i64) -> Vec { + let mut out = vec![C::ZERO; self.amp.len()]; + let x_mask: usize = (0..self.qubit_count) + .filter(|&qubit| x_bits[qubit]) + .map(|qubit| 1usize << (self.qubit_count - 1 - qubit)) + .sum(); + for base in 0..(1 << self.qubit_count) { + let target = base ^ x_mask; + let mut sign_parity = 0i64; + for qubit in 0..self.qubit_count { + if z_bits[qubit] && (base >> (self.qubit_count - 1 - qubit)) & 1 == 1 { + sign_parity ^= 1; + } + } + let coeff = zeta8(2 * phase + 4 * sign_parity); + out[target] += self.amp[base] * coeff; + } + out + } + pub fn apply_pauli(&mut self, x_bits: &[bool], z_bits: &[bool], phase: i64) { + self.amp = self.pauli_applied(x_bits, z_bits, phase); + } + pub fn apply_pauli_exp(&mut self, x_bits: &[bool], z_bits: &[bool], phase: i64) { + let pauli_applied = self.pauli_applied(x_bits, z_bits, phase); + for base in 0..self.amp.len() { + self.amp[base] = (self.amp[base] + pauli_applied[base] * Complex::I) * ROOT_HALF; + } + } + pub fn apply_controlled_pauli( + &mut self, + first_pauli: &(Vec, Vec, i64), + second_pauli: &(Vec, Vec, i64), + ) { + // controlled_pauli(first, second) = (I + first)/2 + (I - first)/2 * second + let first_pauli_applied = self.pauli_applied(&first_pauli.0, &first_pauli.1, first_pauli.2); + let plus: Vec = (0..self.amp.len()) + .map(|index| (self.amp[index] + first_pauli_applied[index]) * 0.5) + .collect(); + let minus = Dense { + qubit_count: self.qubit_count, + amp: (0..self.amp.len()) + .map(|index| (self.amp[index] - first_pauli_applied[index]) * 0.5) + .collect(), + }; + let second_pauli_applied_to_minus = minus.pauli_applied(&second_pauli.0, &second_pauli.1, second_pauli.2); + for index in 0..self.amp.len() { + self.amp[index] = plus[index] + second_pauli_applied_to_minus[index]; + } + } + pub fn project(&mut self, x_bits: &[bool], z_bits: &[bool], phase: i64, outcome: bool) { + let pauli_applied = self.pauli_applied(x_bits, z_bits, phase); + let sign = if outcome { -1.0 } else { 1.0 }; + for index in 0..self.amp.len() { + self.amp[index] = (self.amp[index] + pauli_applied[index] * sign) * 0.5; + } + normalize(&mut self.amp); + } +} + +/// Renormalises an amplitude vector to unit norm. +/// +/// # Panics +/// Panics if the state has (near) zero norm. +pub fn normalize(amp: &mut [C]) { + let norm = amp.iter().map(Complex::norm_sqr).sum::().sqrt(); + assert!(norm > 1e-9, "attempted to normalize a vanishing state"); + let inv = 1.0 / norm; + for amplitude in amp.iter_mut() { + *amplitude *= inv; + } +} + +/// Returns the `2x2` matrix of a single-qubit unitary operation. +/// +/// # Panics +/// Panics if `op` is not a single-qubit operation. +pub fn gate_matrix(op: UnitaryOp) -> [[C; 2]; 2] { + let root_half = ROOT_HALF; + match op { + UnitaryOp::Hadamard => [ + [C::new(root_half, 0.0), C::new(root_half, 0.0)], + [C::new(root_half, 0.0), C::new(-root_half, 0.0)], + ], + UnitaryOp::X => [[C::ZERO, C::new(1.0, 0.0)], [C::new(1.0, 0.0), C::ZERO]], + UnitaryOp::Y => [[C::ZERO, C::new(0.0, -1.0)], [C::new(0.0, 1.0), C::ZERO]], + UnitaryOp::Z => [[C::new(1.0, 0.0), C::ZERO], [C::ZERO, C::new(-1.0, 0.0)]], + UnitaryOp::SqrtZ => [[C::new(1.0, 0.0), C::ZERO], [C::ZERO, C::new(0.0, 1.0)]], + UnitaryOp::SqrtZInv => [[C::new(1.0, 0.0), C::ZERO], [C::ZERO, C::new(0.0, -1.0)]], + UnitaryOp::SqrtX => [ + [zeta8(1) * root_half, zeta8(7) * root_half], + [zeta8(7) * root_half, zeta8(1) * root_half], + ], + UnitaryOp::SqrtXInv => [ + [zeta8(7) * root_half, zeta8(1) * root_half], + [zeta8(1) * root_half, zeta8(7) * root_half], + ], + UnitaryOp::SqrtY => [ + [zeta8(1) * root_half, zeta8(5) * root_half], + [zeta8(1) * root_half, zeta8(1) * root_half], + ], + UnitaryOp::SqrtYInv => [ + [zeta8(7) * root_half, zeta8(7) * root_half], + [zeta8(3) * root_half, zeta8(7) * root_half], + ], + other => panic!("gate_matrix called on multi-qubit op {other:?}"), + } +} + +/// Materialises the dense statevector produced by a phased Clifford unitary +/// acting on `|0...0>`. +pub fn statevector(phased: &PhasedCliffordUnitary) -> Vec { + use binar::BitMatrix; + use binar::matrix::AlignedBitMatrix; + let qubit_count = phased.num_qubits(); + let mut matrix = AlignedBitMatrix::zeros(qubit_count, qubit_count); + for generator in 0..qubit_count { + let image: DensePauli = phased.clifford().image_z(generator); + for qubit in image.x_bits().support() { + matrix.row_mut(generator).assign_index(qubit, true); + } + } + let rank = BitMatrix::from_aligned(matrix).rank(); + let mag = (0.5f64).powf(rank as f64 / 2.0); + let mut out = vec![C::ZERO; 1 << qubit_count]; + for idx in 0..(1usize << qubit_count) { + let mut value = 0usize; + for qubit in 0..qubit_count { + if (idx >> (qubit_count - 1 - qubit)) & 1 == 1 { + value |= 1usize << qubit; + } + } + if let Some(exp) = phased.state_amplitude_phase_exponent_usize(value) { + out[idx] = zeta8(i64::from(exp)) * mag; + } + } + out +} + +/// Decomposes a [`DensePauli`] into `(x_bits, z_bits, xz_phase_exponent)` arrays +/// sized to `qubit_count`. +pub fn pauli_arrays(pauli: &DensePauli, qubit_count: usize) -> (Vec, Vec, i64) { + let mut x_bits = vec![false; qubit_count]; + let mut z_bits = vec![false; qubit_count]; + for qubit in pauli.x_bits().support() { + x_bits[qubit] = true; + } + for qubit in pauli.z_bits().support() { + z_bits[qubit] = true; + } + (x_bits, z_bits, i64::from(pauli.xz_phase_exponent())) +} + +/// Approximate equality of two amplitude vectors, including global phase. +pub fn close(left: &[C], right: &[C]) -> bool { + left.len() == right.len() + && left + .iter() + .zip(right) + .all(|(left_value, right_value)| (left_value - right_value).norm_sqr() < 1e-6) +}