Introduction
Compiling a fault-tolerant quantum algorithm into a concrete circuit requires navigating a large design space of arithmetic decompositions, ancilla-management schedules, and precision trade-offs. Until recently, researchers explored this space by hand, often spending months optimizing a single subroutine and proving its correctness step by step. Automated circuit optimizers [1 , 2 , 3 ] help locally, but they preserve exact unitaries on small subcircuits rather than changing the algorithm or its approximations. AI coding agents [4 ] can automate this search if a fast test harness checks each candidate circuit and scores its resource cost in a closed loop. For general quantum circuits, checking whether a candidate circuit matches a target unitary is QMA-complete [5 ] and intractable on a classical computer above roughly 40 qubits. Fortunately, in many fault-tolerant quantum algorithms, the subroutines that dominate the cost do not use arbitrary quantum interference internally [6 , 7 , 8 , 9 , 10 , 11 , 12 ] . Instead, they use only reversible classical gates, diagonal phase gates, and measurement-based uncomputation. This structure lets a classical simulator efficiently check their approximate correctness on sampled computational basis inputs.
Ref. [13 ] used this method in March 2026 to estimate the cost of breaking 256-bit elliptic curve cryptography on the secp256k1 curve with Shor’s algorithm. To follow responsible disclosure, the authors withheld the explicit circuits. Instead, they released an open-source simulator and test harness with a zero-knowledge proof that their hidden circuit passed the randomized tests. Soon after Ref. [13 ] appeared, the ecdsa.fail challenge took their simulator and test harness and built a public competition for secp256k1 point addition circuits. The harness checked the approximate correctness of each submission in seconds and ranked it by its qubit and Toffoli counts. This fast verification let the organizers accept submissions without trusting the participants. Over 100 participants, most without a quantum computing background, used AI coding agents to further reduce the Toffoli counts by roughly half [14 ] .
In this paper, we apply this workflow to three benchmarks beyond cryptanalysis and use them to seed a Kaggle competition. We make three main contributions:
We release kickmix, an open-source C++20 and Python package in Qualtran [15 ] that builds and simulates circuits with over 10 8 10^{8} 1 0 8
We describe a non-interactive verifier that seeds its test inputs with a salted hash of the submitted circuit. We also bound the distance between the implemented and ideal output states by the circuit’s failure rate on those inputs.
We define three benchmarks and provide a reference implementation for each of them:
Phasing by 1 / x 1/\sqrt{x} 1/ x [9 ] .
Reed-Solomon syndrome decoding for decoded quantum interferometry (DQI) on the optimal polynomial intersection problem [12 , 16 ] .
Quantum lookup tables and Majorana rotations for block encodings of electronic-structure Hamiltonians [11 ] .
Because classical fuzz testing applies only to kickmix circuits, we restrict each benchmark to the kickmix subroutine that dominates the algorithm’s cost. When that subroutine interacts with non-kickmix parts such as arbitrary state preparation, as in Benchmark III, we price those parts with analytical cost formulas. This pricing ensures that the benchmark score still reflects the full algorithm.
Section II defines the kickmix gate set, describes the kickmix software package, and derives the state-distance bound from fuzz testing. Section III describes the background, reference circuit, test distribution, and score for each benchmark. Section IV discusses how circuit testing might change as AI agents improve. The appendices derive the DFTHC block-encoding circuits and the polynomial approximations in the 1 / x 1/\sqrt{x} 1/ x
Kickmix circuits and fuzz testing via Fiat-Shamir heuristic
Optimizing quantum circuits with AI agents requires a fast, reliable way to check whether a candidate circuit works. For general quantum circuits, classical simulation takes exponential time, and deciding whether two circuits implement the same operation is coNP-complete classically and QMA-complete quantumly [5 ] . Ref. [13 ] introduced a restricted family of quantum circuits, called kickmix circuits, built exclusively from reversible classical gates, diagonal phase operations, and measurement-based uncomputation. Because every instruction in this family maps any computational basis state to a single computational basis state up to a scalar phase, a classical simulator can run a kickmix circuit on computational basis inputs in time linear in gate count.
Although full quantum algorithms rely on quantum interference across superpositions, many of the subroutines that dominate fault-tolerant resource estimates (such as integer arithmetic, quantum lookup tables, and syndrome decoders) can be written as kickmix circuits [9 , 12 , 13 , 17 ] . This structure suggests checking candidate compilations by testing them on randomly sampled classical inputs. However, doing so requires care: one must check in each problem context that approximate correctness on sampled basis states guarantees high fidelity on the quantum states used by the algorithm. When this holds, combining classical trajectory simulation with the Fiat-Shamir heuristic [13 , 18 ] turns random-input fuzz testing into a public, non-interactive verification protocol that requires no secret test set. To run this verification loop on subroutines exceeding 10 8 10^8 1 0 8 kickmix as an open-source C++20 and Python package inside Qualtran [15 ] . On a single CPU core, kickmix builds circuits at 2 × 10 7 2 \times10^7 2 × 1 0 7 6.0 × 10 10 6.0 \times10^{10} 6.0 × 1 0 10
The kickmix gate set
When a kickmix circuit runs on an n n n ∣ x ⟩ |x\rangle ∣ x ⟩
Reversible classical gates (X, CX, CCX, and SWAP): Permute computational basis states by updating the stored qubit bits directly.
Diagonal phase gates (Z, CZ, CCZ, Z_POW, and NEG): Leave the qubit bits unchanged and update the tracked output phase whenever all control qubits hold ∣ 1 ⟩ |1\rangle ∣1 ⟩
Measurement and reset gates (HMR and R): Measure a qubit in the X basis into a classical bit and reset it to ∣ 0 ⟩ |0\rangle ∣0 ⟩ ∣ 0 ⟩ |0\rangle ∣0 ⟩
Classical bit operations (BIT_STORE0, BIT_STORE1, and BIT_INVERT): Set classical register bits to 0 or 1, or invert them in place.
Every instruction above can be controlled per gate (if b) or across blocks via a condition stack (PUSH_CONDITION and POP_CONDITION), running only when all active classical conditions hold 1 1 1
The HMR instruction (pronounced “hammer”) applies a Hadamard gate, measures the target qubit q t q_t q t Z Z Z q t q_t q t ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ q t q_t q t f ( x ) f(x) f ( x ) X X X b ∈ { 0 , 1 } b \in\{0,1\} b ∈ { 0 , 1 }
∣ x ⟩ ∣ f ( x ) ⟩ → H M R ( − 1 ) b ⋅ f ( x ) ∣ x ⟩ ∣ 0 ⟩ . (1) \lvert x\rangle\lvert f(x)\rangle\xrightarrow{\mathrm{HMR}} (-1)^{b\cdot f(x)}\lvert x\rangle\lvert0\rangle. \tag*{(1)} ∣ x ⟩ ∣ f ( x )⟩ HMR ( − 1 ) b ⋅ f ( x ) ∣ x ⟩ ∣ 0 ⟩ . (1) When b = 1 b=1 b = 1 ( − 1 ) f ( x ) (-1)^{f(x)} ( − 1 ) f ( x ) X X X ∣ x 1 ∧ x 2 ⟩ \lvert x_1 \land x_2\rangle ∣ x 1 ∧ x 2 ⟩ b = 1 b=1 b = 1 [17 ] .
However, HMR only works when the resulting phase error can be corrected cheaply. If a circuit applies HMR to a junk qubit that remains entangled with the system without fixing its phase kickback, classical trajectory simulation fails loudly. The simulator samples a uniformly random bit b ∈ { 0 , 1 } b \in\{0,1\} b ∈ { 0 , 1 } q t ⋅ b = 1 q_t \cdot b=1 q t ⋅ b = 1 q t q_t q t ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ 1 / 2 1/2 1/2 q t = 1 q_t=1 q t = 1 ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ q t = 1 q_t=1 q t = 1
The kickmix software package in Qualtran
The kickmix library is implemented in C++20 with Python bindings inside Qualtran [15 ] . Inspired by the bit-packed stabilizer engine of Stim [19 ] , kickmix is built around two design principles that let humans and AI agents iterate on 10 8 10^8 1 0 8
First, performance drives every layer of circuit construction, serialization, and simulation. Standard Python circuit libraries allocate a separate heap object for every gate and simulate one input at a time, which becomes prohibitively slow at 10 8 10^8 1 0 8 kickmix stores a circuit using a struct-of-arrays representation grouped by operand count: one byte array of opcodes (OpType) alongside separate arrays for 1-, 2-, and 3-qubit arguments, classical target and condition bits, and rotation angles. Each 1-byte opcode embeds bit flags indicating which operand arrays it consumes, letting the simulator step sequential pointers through only the active arrays without padding every gate to the worst-case size. Saving and loading a circuit in binary .km format transfers each array between memory and disk as a contiguous binary block, avoiding per-gate text formatting and parsing. During simulation, kickmix bit-packs 64 to 512 independent trajectories across CPU SIMD registers so each instruction updates an entire batch of shots in parallel. The one exception is Z_POW: because each trajectory maintains its own 128-bit fixed-precision phase accumulator rather than a single sign bit, a Z_POW instruction cannot update the whole batch with a single bitwise SIMD operation. Instead, it scans the active shot mask across all trajectories in the batch and adds the rotation angle to each active shot sequentially.
Second, expressiveness guides the Python builder API so complex end-to-end algorithms remain compact and readable. To eliminate the repetitive boilerplate that arises when writing multi-stage fault-tolerant algorithms with measurement-based uncomputation, kickmix provides four abstractions:
Automatic downgrading: Operations accept qubits, classical bits, and compile-time constants interchangeably, automatically replacing classically controlled gates with lower-cost instructions, dropping no-ops on False controls, and supporting block-scoped classical conditions.
NumPy-style broadcasting: Most primitive gates, measurements, and classical bit operations accept both single wires and register arrays (km.array), broadcasting scalars against registers or operating elementwise across registers, so a single builder call replaces an explicit loop over individual qubits.
Dynamic workspace allocation: The builder provides methods to allocate and free temporary quantum and classical registers within a configurable qubit budget. Subroutines can use these methods to request and release scratch space locally, and freed qubits are handed back out to later allocations.
Library of common subroutines: The builder provides a library of commonly used fault-tolerant subroutines, including temporary-AND gates, QROM table lookups, adders, comparators, finite-field arithmetic, and context-managed unary iteration, with paired compute and measurement-based uncompute variants.
kickmix also includes inspection utilities that generate circuit diagrams and interactive SVG flame charts showing where Toffoli gates and qubits are spent across nested subroutines.
Fuzz testing via the Fiat-Shamir heuristic
To verify a candidate kickmix circuit, we simulate it on computational basis inputs using km.Simulator and compare the resulting bitstrings and phases against a classical reference. When the input domain is small, we can test every basis state exhaustively to certify exact correctness on all inputs. When the input domain is too large to enumerate, fuzz testing on a random sample of inputs instead certifies approximate correctness by bounding the probability ε \varepsilon ε ε \varepsilon ε [7 , 13 ] .
A bound on the per-input failure rate ε \varepsilon ε r ∈ { 0 , 1 } m r \in\{0,1\}^{m} r ∈ { 0 , 1 } m m m m V r V_r V r b b b ( ∣ 0 ⟩ ⟨ 0 ∣ + ( − 1 ) b ∣ 0 ⟩ ⟨ 1 ∣ ) / 2 (|0\rangle\langle0|+(-1)^b|0\rangle\langle1|)/\sqrt{2} ( ∣0 ⟩ ⟨ 0∣ + ( − 1 ) b ∣0 ⟩ ⟨ 1∣ ) / 2 2 − m / 2 V r 2^{-m/2}V_r 2 − m /2 V r r r r 2 − m 2^{-m} 2 − m ∣ ψ ⟩ = ∑ x α x ∣ x ⟩ |\psi\rangle=\sum_x\alpha_x|x\rangle ∣ ψ ⟩ = ∑ x α x ∣ x ⟩
W ∣ ψ ⟩ = 2 − m / 2 ∑ r V r ∣ ψ ⟩ ⊗ ∣ r ⟩ . (2) W|\psi\rangle=2^{-m/2}\sum_r V_r|\psi\rangle\otimes|r\rangle. \tag*{(2)} W ∣ ψ ⟩ = 2 − m /2 r ∑ V r ∣ ψ ⟩ ⊗ ∣ r ⟩ . (2) The ideal circuit applies the target unitary U U U ∣ + ⟩ ⊗ m |+\rangle^{\otimes m} ∣ + ⟩ ⊗ m B r = { x : V r ∣ x ⟩ ≠ U ∣ x ⟩ } B_r=\{x:V_r|x\rangle\ne U|x\rangle\} B r = { x : V r ∣ x ⟩ = U ∣ x ⟩} r r r ε r = ∑ x ∈ B r ∣ α x ∣ 2 \varepsilon_r=\sum_{x\in B_r}|\alpha_x|^2 ε r = ∑ x ∈ B r ∣ α x ∣ 2 ε = E r [ ε r ] \varepsilon=\mathbb{E}_r[\varepsilon_r] ε = E r [ ε r ] x ∼ ∣ α x ∣ 2 x\sim|\alpha_x|^2 x ∼ ∣ α x ∣ 2 r r r ∣ 1 ⟩ |1\rangle ∣1 ⟩ b = 1 b=1 b = 1 b = 0 b=0 b = 0 ε \varepsilon ε b = 0 b=0 b = 0
We assume that each V r V_r V r ∣ ψ ⟩ |\psi\rangle ∣ ψ ⟩ V r V_r V r ∣ ψ ⟩ |\psi\rangle ∣ ψ ⟩ U U U B r B_r B r ∥ ( V r − U ) ∑ x ∈ B r α x ∣ x ⟩ ∥ ≤ 2 ε r \|(V_r-U)\sum_{x\in B_r}\alpha_x|x\rangle\|\le2\sqrt{\varepsilon_r} ∥ ( V r − U ) ∑ x ∈ B r α x ∣ x ⟩ ∥ ≤ 2 ε r r r r
∥ W ∣ ψ ⟩ − U ∣ ψ ⟩ ⊗ ∣ + ⟩ ⊗ m ∥ 2 = E r [ ∥ ( V r − U ) ∑ x ∈ B r α x ∣ x ⟩ ∥ 2 ] ≤ 4 ε . (3) \|W|\psi\rangle-U|\psi\rangle\otimes|+\rangle^{\otimes m}\|^2=\mathbb{E}_r\left[\left\|(V_r-U)\sum_{x\in B_r}\alpha_x|x\rangle\right\|^2\right]\le4\varepsilon. \tag*{(3)} ∥ W ∣ ψ ⟩ − U ∣ ψ ⟩ ⊗ ∣ + ⟩ ⊗ m ∥ 2 = E r ( V r − U ) x ∈ B r ∑ α x ∣ x ⟩ 2 ≤ 4 ε . (3) Two unit vectors at Euclidean distance d d d 1 − d 2 / 2 1-d^2/2 1 − d 2 /2 ( 1 − 2 ε ) 2 (1-2\varepsilon)^2 ( 1 − 2 ε ) 2 U ∣ ψ ⟩ ⊗ ∣ + ⟩ ⊗ m U|\psi\rangle\otimes|+\rangle^{\otimes m} U ∣ ψ ⟩ ⊗ ∣ + ⟩ ⊗ m 1 − ( 1 − 2 ε ) 2 ≤ 2 ε \sqrt{1-(1-2\varepsilon)^2}\le2\sqrt{\varepsilon} 1 − ( 1 − 2 ε ) 2 ≤ 2 ε U ∣ ψ ⟩ U|\psi\rangle U ∣ ψ ⟩ ∣ α x ∣ 2 |\alpha_x|^2 ∣ α x ∣ 2 1
To prevent a submission from hard-coding outputs on a fixed test set, our verifier generates challenge inputs via the Fiat-Shamir heuristic [13 , 18 ] . Given a serialized circuit file C C C C C C C C C N N N
On each shot i ∈ { 1 , … , N } i\in\{1,\ldots,N\} i ∈ { 1 , … , N }
Input and output registers hold their expected classical bitstrings.
Every scratch qubit returns cleanly to ∣ 0 ⟩ |0\rangle ∣0 ⟩
The accumulated output phase matches the target phase within the benchmark tolerance.
If a shot violates any of these conditions, the verifier records a failure. Suppose a circuit fails on k k k N N N ε ^ = k / N \hat{\varepsilon}=k/N ε ^ = k / N ε \varepsilon ε ε + > ε ^ \varepsilon^{+}>\hat{\varepsilon} ε + > ε ^ k k k N N N
Pr ( at most k failures in N shots ∣ ε ≥ ε + ) ≤ ∑ j = 0 k ( N j ) ( ε + ) j ( 1 − ε + ) N − j . (4) \Pr(\text{at most } k \text{ failures in } N \text{ shots} \mid\varepsilon\geq\varepsilon^{+}) \leq\sum_{j=0}^{k} \binom{N}{j}(\varepsilon^{+})^{j}(1-\varepsilon^{+})^{N-j}. \tag*{(4)} Pr ( at most k failures in N shots ∣ ε ≥ ε + ) ≤ j = 0 ∑ k ( j N ) ( ε + ) j ( 1 − ε + ) N − j . (4) Setting the right-hand side of Eq. (4) equal to a target false-pass probability δ = 2 − 128 \delta=2^{-128} δ = 2 − 128 ε + \varepsilon^{+} ε + ε \varepsilon ε 1 − 2 − 128 1-2^{-128} 1 − 2 − 128 N N N k = 0 k=0 k = 0 j = 0 j=0 j = 0 ( 1 − ε + ) N = δ (1-\varepsilon^{+})^{N}=\delta ( 1 − ε + ) N = δ ε + ≈ ln ( 1 / δ ) / N ≈ 88.7 / N \varepsilon^{+}\approx\ln(1/\delta)/N\approx88.7/N ε + ≈ ln ( 1/ δ ) / N ≈ 88.7/ N N = 10,240 N=10{,}240 N = 10 , 240 ε + = 1 % \varepsilon^{+}=1\% ε + = 1% 0.99 10240 < 2 − 128 0.99^{10240}<2^{-128} 0.9 9 10240 < 2 − 128 N = 102,400 N=102{,}400 N = 102 , 400 ε + < 0.1 % \varepsilon^{+}<0.1\% ε + < 0.1% 1 − 2 − 128 1-2^{-128} 1 − 2 − 128 k = 0 k=0 k = 0 N N N k ≥ 0 k\geq0 k ≥ 0 ε + \varepsilon^{+} ε + 2
Kaggle competition
We are launching a public competition on Kaggle, an online platform for machine learning and algorithmic challenges, featuring three benchmarks from fusion simulation, decoded quantum interferometry, and quantum chemistry. Each benchmark isolates the kickmix subroutine that dominates the fault-tolerant cost of the full algorithm at a fixed problem size. In our own runs with AI coding agents, we have found significant improvements over the previously published baselines, and we expect further improvements remain to be discovered. For each benchmark, we commit to our best internal run via its SHA-256 hash2
Benchmark I: Inverse square roots for Coulomb simulations with relevance to fusion reactor design
Problem introduction and importance
Quantum simulation of matter is one of the primary applications for quantum computers. This simulation capability is even more important where experimental verification is challenging. One such example we recently explored was the calculation of stopping power relevant to warm-dense-matter simulation [9 ] . The calculation of a stopping power is critical to the understanding and design of inertial confinement fusion devices. The stopping power calculation itself is representative of a broader class of dynamical simulations containing core quantum circuit primitives. Thus our focus here is on one of the highest cost algorithmic building blocks. The primitive we focus on in this section is the primary bottleneck for the Trotterized time-evolution for a system represented in a first-quantized representation [9 ] and involves the coherent calculation of the inverse distance between two particles.
Our goal will be to implement the dynamics of η \eta η
H = T + U + V (5) H = T + U + V \tag*{(5)} H = T + U + V (5) which involves three terms: 1) T T T U U U V V V U U U V V V e − i δ T e^{-i\delta T} e − i δ T e − i δ ( U + V ) e^{-i\delta(U+V)} e − i δ ( U + V ) δ \delta δ
U = ∑ i = 1 η ∑ ℓ = 1 L Z ℓ ∣ r i − R ℓ ∣ = ∑ i = 1 η ∑ ℓ = 1 L Z ℓ ( x i − R x , ℓ ) 2 + ( y i − R y , ℓ ) 2 + ( z i − R z , ℓ ) 2 (6) U = \sum_{i=1}^{\eta}\sum_{\ell=1}^{L}\frac{Z_{\ell}}{|r_i-R_{\ell}|} = \sum_{i=1}^{\eta}\sum_{\ell=1}^{L}\frac{Z_{\ell}}{\sqrt{(x_i-R_{x,\ell})^2+(y_i-R_{y,\ell})^2+(z_i-R_{z,\ell})^2}} \tag*{(6)} U = i = 1 ∑ η ℓ = 1 ∑ L ∣ r i − R ℓ ∣ Z ℓ = i = 1 ∑ η ℓ = 1 ∑ L ( x i − R x , ℓ ) 2 + ( y i − R y , ℓ ) 2 + ( z i − R z , ℓ ) 2 Z ℓ (6) V = ∑ i , j = 1 i ≠ j η 1 ∣ r i − r j ∣ = ∑ i , j = 1 i ≠ j η 1 ( x i − x j ) 2 + ( y i − y j ) 2 + ( z i − z j ) 2 (7) V = \sum_{\substack{i,j=1\\i\ne j}}^{\eta}\frac{1}{|r_i-r_j|} = \sum_{\substack{i,j=1\\i\ne j}}^{\eta}\frac{1}{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2}} \tag*{(7)} V = i , j = 1 i = j ∑ η ∣ r i − r j ∣ 1 = i , j = 1 i = j ∑ η ( x i − x j ) 2 + ( y i − y j ) 2 + ( z i − z j ) 2 1 (7) where we have expanded distance in the Euclidean directions to indicate the main computation needed is the inverse square root of the sum of squared differences in each direction. Our focus will be the electron-electron computation as this needs to be performed ( η 2 ) \binom{\eta}{2} ( 2 η ) η \eta η n p n_p n p
∣ ψ ⟩ = ∑ x , y , z = 0 2 n p − 1 c ( x , y , z ) ∣ x ⟩ ∣ y ⟩ ∣ z ⟩ = ∑ r c ( r ) ∣ r ⟩ . (8) |\psi\rangle= \sum_{x,y,z=0}^{2^{n_p}-1} c(x,y,z)|x\rangle|y\rangle|z\rangle= \sum_r c(r)|r\rangle. \tag*{(8)} ∣ ψ ⟩ = x , y , z = 0 ∑ 2 n p − 1 c ( x , y , z ) ∣ x ⟩ ∣ y ⟩ ∣ z ⟩ = r ∑ c ( r ) ∣ r ⟩ . (8) The quantum algorithm involves pairwise computation of V ( r i , r j ) V(r_i,r_j) V ( r i , r j ) e − i δ V ( r i , r j ) e^{-i\delta V(r_i,r_j)} e − i δ V ( r i , r j ) V ( r i , r j ) V(r_i,r_j) V ( r i , r j ) [9 ] performing the computation pairwise instead of for all pairs and summing saves ancilla qubits. The phasing algorithm looks like the following (reproduced from Ref. [20 ] , Eqs. (15) –(17) ):
∑ r 1 , … , r η c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ 000 … ⟩ (9) \sum_{r_1,\ldots,r_\eta} c(r_1,\ldots,r_\eta)|r_1,\ldots,r_\eta\rangle|000\ldots\rangle \tag*{(9)} r 1 , … , r η ∑ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣000 … ⟩ (9) → ∑ r 1 , … , r η c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ V ( r i , r j ) ⟩ (10) \rightarrow\sum_{r_1,\ldots,r_\eta} c(r_1,\ldots,r_\eta)|r_1,\ldots,r_\eta\rangle|V(r_i,r_j)\rangle \tag*{(10)} → r 1 , … , r η ∑ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ V ( r i , r j )⟩ (10) → ∑ r 1 , … , r η e − i V ( r i , r j ) δ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ V ( r i , r j ) ⟩ (11) \rightarrow\sum_{r_1,\ldots,r_\eta}e^{-iV(r_i,r_j)\delta}c(r_1,\ldots,r_\eta)|r_1,\ldots,r_\eta\rangle|V(r_i,r_j)\rangle \tag*{(11)} → r 1 , … , r η ∑ e − iV ( r i , r j ) δ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ V ( r i , r j )⟩ (11) → ∑ r 1 , … , r η e − i V ( r i , r j ) δ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣ 000 … ⟩ (12) \rightarrow\sum_{r_1,\ldots,r_\eta}e^{-iV(r_i,r_j)\delta}c(r_1,\ldots,r_\eta)|r_1,\ldots,r_\eta\rangle|000\ldots\rangle \tag*{(12)} → r 1 , … , r η ∑ e − iV ( r i , r j ) δ c ( r 1 , … , r η ) ∣ r 1 , … , r η ⟩ ∣000 … ⟩ (12) The strategy adopted in Refs. [9 , 21 ] follows a similar trick that was used in the video game “Quake 3” employing a single step of Newton-Raphson after construction of a high enough quality initial guess. The C code [22 ] for this is included below
float Q_rsqrt( float number )
{
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = * ( long * ) &y; // evil floating point bit level hacking
i = 0x5f3759df - ( i >> 1 ); // what the f***?
y = * ( float * ) &i;
y = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
#ifndef Q3_VM
#ifdef __linux__
assert ( !isnan(y) ); // bk010122 - FPE?
#endif
#endif
return y;
}which produces an approximation that has a maximum relative error of 0.00175228 over the domain of 32 bit floats (IEEE 754-1985). We work with fixed point representation on quantum computers so we will design an algorithm that works approximately the same but with a fixed point.
Due to the accuracy requirements of the total time evolution Ref. [9 ] determined that 24-bits of relative error would be needed for the computation in the inverse square root. Starting from an input register ∣ x ⟩ \lvert x\rangle ∣ x ⟩ x x x [15 ] . Our primary task for an AI agent to optimize will be the coherent arithmetic and phasing costs for computing steps 1–3 while maintaining the 24-bit relative error precision bound. Our baseline implementation is merely a starting point for optimization as there are a number of choices one can make in the computation of the inverse square root. For example, not using a cubic polynomial over a small rescaled range but a simpler polynomial over a larger range may reduce cost further. Also the accuracy of coefficients output by QROM can be further optimized. One could imagine going a completely different route and simply finding a good polynomial approximation over the input domain. The number of tradeoffs and directions to explore (as in many compilation tasks) makes circuit primitive optimization of this type well posed for AI agent exploration. Restricting to kickmix circuits allows efficient exploration of these optimizations by many AI agents and will likely lead to large space-time cost reductions.
Reference implementation
Calculation of the inverse square-root was sketched in Ref. [9 ] in only enough detail to justify the constant factor analysis. Here we describe the full implementation which models the classical fixed point method using a single step of Newton-Raphson and provide the reference implementation of the paper. Computing w = ( x i − x j ) 2 + ( y i − y j ) 2 + ( z i − z j ) 2 w=(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2 w = ( x i − x j ) 2 + ( y i − y j ) 2 + ( z i − z j ) 2 n p n_p n p 3 n p 2 − n p − 1 3n_p^2-n_p-1 3 n p 2 − n p − 1 [23 ] ) and lies upstream of the benchmark subroutine, which takes ∣ x ⟩ = ∣ w ⟩ \lvert x\rangle=\lvert w\rangle ∣ x ⟩ = ∣ w ⟩ ( 2 n p + 2 ) (2n_p+2) ( 2 n p + 2 ) w w w 1 / w n o r m 1/\sqrt{w_{\mathrm{norm}}} 1/ w norm n p = 6 n_p=6 n p = 6 w w w [ 1 , 11532 ] [1,11532] [ 1 , 11532 ] n p = 10 n_p=10 n p = 10 η = 218 \eta=218 η = 218 [ 1 , 3 ( 2 10 − 2 ) 2 ] = [ 1 , 3 , 133 , 452 ] [1,3(2^{10}-2)^2]=[1,3,133,452] [ 1 , 3 ( 2 10 − 2 ) 2 ] = [ 1 , 3 , 133 , 452 ] w w w 1 / w 1/\sqrt{w} 1/ w w = 1 ↦ 1 / w = 1 w=1\mapsto1/\sqrt{w}=1 w = 1 ↦ 1/ w = 1 w = 11532 ↦ 1 / w ≈ 0.0093 w=11532\mapsto1/\sqrt{w}\approx0.0093 w = 11532 ↦ 1/ w ≈ 0.0093 ∼ 7 \sim7 ∼ 7 ∼ 11 \sim11 ∼ 11 n p = 10 n_p=10 n p = 10 n p = 6 n_p=6 n p = 6 24 + 7 = 31 24+7=31 24 + 7 = 31 n p = 10 n_p=10 n p = 10
First, Refs. [9 , 21 ] avoid the high multiplication cost by normalizing w w w [ 1 , 4 ) [1,4) [ 1 , 4 ) w w w n p = 6 n_p=6 n p = 6
w = 3 ↦ 00 00 00 00 00 00 11 (12 leading zeros) w = 17 ↦ 00 00 00 00 01 00 01 (8 leading zeros) w = 4096 ↦ 01 00 00 00 00 00 00 (1 leading zero) w = 11532 ↦ 10 11 01 00 00 11 00 (0 leading zeros) (13) \begin{aligned}
w=3 &\quad\mapsto00\ 00\ 00\ 00\ 00\ 00\ 11\ \text{(12 leading zeros)}\\
w=17 &\quad\mapsto00\ 00\ 00\ 00\ 01\ 00\ 01\ \text{(8 leading zeros)}\\
w=4096 &\quad\mapsto01\ 00\ 00\ 00\ 00\ 00\ 00\ \text{(1 leading zero)}\\
w=11532 &\quad\mapsto10\ 11\ 01\ 00\ 00\ 11\ 00\ \text{(0 leading zeros)} \tag*{(13)}
\end{aligned} w = 3 w = 17 w = 4096 w = 11532 ↦ 00 00 00 00 00 00 11 (12 leading zeros) ↦ 00 00 00 00 01 00 01 (8 leading zeros) ↦ 01 00 00 00 00 00 00 (1 leading zero) ↦ 10 11 01 00 00 11 00 (0 leading zeros) (13) all these numbers have different significant figures and thus padding all the values to the appropriate significant figures for 2 − 24 2^{-24} 2 − 24 w w w
w n o r m = w ⋅ 2 2 s (14) w_{\mathrm{norm}}=w\cdot2^{2s} \tag*{(14)} w norm = w ⋅ 2 2 s (14) where s s s w ↦ w n o r m w\mapsto w_{\mathrm{norm}} w ↦ w norm w n o r m ∈ [ 1 , 4 ) w_{\mathrm{norm}}\in[1,4) w norm ∈ [ 1 , 4 ) 1 / w n o r m ∈ ( 0.5 , 1 ] 1/\sqrt{w_{\mathrm{norm}}}\in(0.5,1] 1/ w norm ∈ ( 0.5 , 1 ] w n o r m w_{\mathrm{norm}} w norm s s s s s s x x x x n o r m x_{\mathrm{norm}} x norm x x x 2 n p + 2 2n_p+2 2 n p + 2 n p = 6 n_p=6 n p = 6 n p = 10 n_p=10 n p = 10 s s s n p n_p n p ∣ s ⟩ \lvert s\rangle ∣ s ⟩ x x x s s s x x x
After shifting, w n o r m w_{\mathrm{norm}} w norm n p = 6 n_p = 6 n p = 6 w n o r m ∈ [ 4096 , 16383 ] w_{\mathrm{norm}} \in[4096, 16383] w norm ∈ [ 4096 , 16383 ] 1 / w n o r m 1/\sqrt{w_{\mathrm{norm}}} 1/ w norm [ 1 , 4 ) [1,4) [ 1 , 4 ) [ 1 , 2 ] [1,2] [ 1 , 2 ] x n o r m ∈ [ 2 , 4 ) x_{\mathrm{norm}} \in[2,4) x norm ∈ [ 2 , 4 )
We will split the interval [ 1 , 2 ] [1,2] [ 1 , 2 ]
1 x ≈ P a ( x ) = a 0 − a 1 ( x − 1 ) + a 2 ( x − 1 ) 2 − a 3 ( x − 1 ) 3 (15) \frac{1}{\sqrt{x}} \approx P_a(x) = a_0 - a_1(x-1) + a_2(x-1)^2 - a_3(x-1)^3 \tag*{(15)} x 1 ≈ P a ( x ) = a 0 − a 1 ( x − 1 ) + a 2 ( x − 1 ) 2 − a 3 ( x − 1 ) 3 (15) for x ∈ [ 1 , 3 / 2 ] x \in[1,3/2] x ∈ [ 1 , 3/2 ]
1 x ≈ P b ( x ) = b 0 − b 1 ( x − 1.5 ) + b 2 ( x − 1.5 ) 2 − b 3 ( x − 1.5 ) 3 (16) \frac{1}{\sqrt{x}} \approx P_b(x) = b_0 - b_1(x-1.5) + b_2(x-1.5)^2 - b_3(x-1.5)^3 \tag*{(16)} x 1 ≈ P b ( x ) = b 0 − b 1 ( x − 1.5 ) + b 2 ( x − 1.5 ) 2 − b 3 ( x − 1.5 ) 3 (16) for x ∈ [ 3 / 2 , 2 ] x \in[3/2,2] x ∈ [ 3/2 , 2 ] [9 ] where
a 0 = 0.99994132489119882162 , (17) a_0 = 0.99994132489119882162, \tag*{(17)} a 0 = 0.99994132489119882162 , (17) a 1 = 0.49609891915903542303 , (18) a_1 = 0.49609891915903542303, \tag*{(18)} a 1 = 0.49609891915903542303 , (18) a 2 = 0.33261112772430493331 , (19) a_2 = 0.33261112772430493331, \tag*{(19)} a 2 = 0.33261112772430493331 , (19) a 3 = 0.14876762006038398086. (20) a_3 = 0.14876762006038398086. \tag*{(20)} a 3 = 0.14876762006038398086. (20) and
b 0 = 0.81648515205385221995 , (21) b_0 = 0.81648515205385221995, \tag*{(21)} b 0 = 0.81648515205385221995 , (21) b 1 = 0.27136515484240234115 , (22) b_1 = 0.27136515484240234115, \tag*{(22)} b 1 = 0.27136515484240234115 , (22) b 2 = 0.12756148214815175348 , (23) b_2 = 0.12756148214815175348, \tag*{(23)} b 2 = 0.12756148214815175348 , (23) b 3 = 0.044753028579153842218. (24) b_3 = 0.044753028579153842218. \tag*{(24)} b 3 = 0.044753028579153842218. (24) In Figure 1 we visualize what this initial guess looks like over the unscaled interval [ 1 , 2 ] [1,2] [ 1 , 2 ] 2 − 24 2^{-24} 2 − 24 [ 2 , 4 ) [2,4) [ 2 , 4 )
Figure 1. Polynomial approximation on [ 1 , 2 ] [1,2] [ 1 , 2 ] [ 1 , 1.5 ] [1,1.5] [ 1 , 1.5 ] [ 1.5 , 2 ] [1.5,2] [ 1.5 , 2 ] { a } \{a\} { a } { b } \{b\} { b } [9 ] .
A full description of the algorithmic primitives for QROM output of the multiplication coefficients, how to rescale the coefficients for the [ 2 , 4 ) [2,4) [ 2 , 4 )
Next, we can see what happens to this w n o r m w_{\mathrm{norm}} w norm w n o r m w_{\mathrm{norm}} w norm
y 0 ≈ 1 / w n o r m (25) y_0 \approx1/\sqrt{w_{\mathrm{norm}}} \tag*{(25)} y 0 ≈ 1/ w norm (25) which is related to our true inverse square root value
y 0 ≈ ( w ⋅ 2 2 s ) − 1 / 2 = w − 1 / 2 ⋅ 2 − s (26) y_0 \approx\left(w \cdot2^{2s}\right)^{-1/2}=w^{-1/2}\cdot2^{-s} \tag*{(26)} y 0 ≈ ( w ⋅ 2 2 s ) − 1/2 = w − 1/2 ⋅ 2 − s (26) if we want to get the “true” initial guess we would need to scale y 0 y_0 y 0 2 s 2^s 2 s
y i n i t = y 0 ⋅ 2 s ≈ 1 / w (27) y_{\mathrm{init}}=y_0\cdot2^s\approx1/\sqrt{w} \tag*{(27)} y init = y 0 ⋅ 2 s ≈ 1/ w (27) We can now do substitutions into our Newton-Raphson step
y 1 = 1 2 y i n i t ( 3 − w ⋅ y i n i t 2 ) (28) y_1=\frac{1}{2}y_{\mathrm{init}}\left(3-w\cdot y_{\mathrm{init}}^2\right) \tag*{(28)} y 1 = 2 1 y init ( 3 − w ⋅ y init 2 ) (28) = 1 2 ( y 0 ⋅ 2 s ) ( 3 − w ⋅ ( y 0 ⋅ 2 s ) 2 ) (29) =\frac{1}{2}\left(y_0\cdot2^s\right)\left(3-w\cdot\left(y_0\cdot2^s\right)^2\right) \tag*{(29)} = 2 1 ( y 0 ⋅ 2 s ) ( 3 − w ⋅ ( y 0 ⋅ 2 s ) 2 ) (29) = 1 2 ( y 0 ⋅ 2 s ) ( 3 − w ⋅ 2 2 s ⋅ y 0 2 ) (30) =\frac{1}{2}\left(y_0\cdot2^s\right)\left(3-w\cdot2^{2s}\cdot y_0^2\right) \tag*{(30)} = 2 1 ( y 0 ⋅ 2 s ) ( 3 − w ⋅ 2 2 s ⋅ y 0 2 ) (30) = 1 2 ( y 0 ⋅ 2 s ) ( 3 − w n o r m ⋅ y 0 2 ) (31) =\frac{1}{2}\left(y_0\cdot2^s\right)\left(3-w_{\mathrm{norm}}\cdot y_0^2\right) \tag*{(31)} = 2 1 ( y 0 ⋅ 2 s ) ( 3 − w norm ⋅ y 0 2 ) (31) which shows we do not need to know the unnormalized w w w w n o r m w_{\mathrm{norm}} w norm s s s w n o r m ⋅ y 0 2 w_{\mathrm{norm}}\cdot y_0^2 w norm ⋅ y 0 2 3 − 1 3-1 3 − 1 w w w w n o r m w_{\mathrm{norm}} w norm y 0 y_0 y 0 p = 1 2 y 0 ( 3 − w n o r m ⋅ y 0 2 ) p=\frac{1}{2}y_0\left(3-w_{\mathrm{norm}}\cdot y_0^2\right) p = 2 1 y 0 ( 3 − w norm ⋅ y 0 2 ) y 0 y_0 y 0 y 1 = p ≪ s y_1=p\ll s y 1 = p ≪ s n p n_p n p 30 n p 30n_p 30 n p y 1 y_1 y 1 ∣ x ⟩ ∣ 0 ⟩ \lvert x\rangle\lvert0\rangle ∣ x ⟩ ∣ 0 ⟩
∣ x ⟩ → shift by 2 s x n o r m → QROM + Cubic Poly y 0 → NR step p = 1 2 y 0 ( 3 − x n o r m y 0 2 ) → ≪ s y 1 = p ⋅ 2 s → Z_POW + uncomp. e 2 π i / x ∣ x ⟩ \lvert x\rangle\xrightarrow{\text{shift by }2^s} x_{\mathrm{norm}}
\xrightarrow{\text{QROM + Cubic Poly}} y_0
\xrightarrow{\text{NR step}} p=\frac{1}{2}y_0\left(3-x_{\mathrm{norm}}y_0^2\right)
\xrightarrow{\ll s} y_1=p\cdot2^s
\xrightarrow{\text{Z\_POW + uncomp.}} e^{2\pi i/\sqrt{x}}\lvert x\rangle ∣ x ⟩ shift by 2 s x norm QROM + Cubic Poly y 0 NR step p = 2 1 y 0 ( 3 − x norm y 0 2 ) ≪ s y 1 = p ⋅ 2 s Z_POW + uncomp. e 2 π i / x ∣ x ⟩ For the benchmark size n p = 10 n_p=10 n p = 10 x x x Q = 260 Q=260 Q = 260 T T T C = 7 , 988 C=7,988 C = 7 , 988 Q × C = 2.08 × 10 6 Q\times C=2.08\times10^6 Q × C = 2.08 × 1 0 6 2.84 × 10 − 8 = 2 − 25.07 2.84\times10^{-8}=2^{-25.07} 2.84 × 1 0 − 8 = 2 − 25.07 x = 2 , 588 , 411 x=2,588,411 x = 2 , 588 , 411
Verification strategy
A submission is a single .kmb circuit file implementing ∣ x ⟩ ∣ 0 ⟩ ↦ exp ( 2 π i β / x ) ∣ x ⟩ ∣ 0 ⟩ \lvert x\rangle\lvert0\rangle\mapsto\exp(2\pi i\beta/\sqrt{x})\lvert x\rangle\lvert0\rangle ∣ x ⟩ ∣ 0 ⟩ ↦ exp ( 2 π i β / x ) ∣ x ⟩ ∣ 0 ⟩ x ∈ [ 1 , max x ] x\in[1,\max x] x ∈ [ 1 , max x ] n p = 10 n_p=10 n p = 10 β = 1 \beta=1 β = 1 2 n p + 2 = 22 2n_p+2=22 2 n p + 2 = 22 max x = 3 ( 2 10 − 2 ) 2 = 3 , 133 , 452 \max x=3(2^{10}-2)^2=3,133,452 max x = 3 ( 2 10 − 2 ) 2 = 3 , 133 , 452 ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ ∣ x ⟩ \lvert x\rangle ∣ x ⟩ x ∈ [ 1 , max x ] x\in[1,\max x] x ∈ [ 1 , max x ] θ ( x ) \theta(x) θ ( x ) 2 − 64 2^{-64} 2 − 64 [ 2 65 β / x ] = isqrt ( [ 2 130 β 2 / x ] ) [2^{65}\beta/\sqrt{x}]=\operatorname{isqrt}([2^{130}\beta^2/x]) [ 2 65 β / x ] = isqrt ([ 2 130 β 2 / x ]) k = 0 k=0 k = 0 tol ( x ) = 2 − 24 β / x \operatorname{tol}(x)=2^{-24}\beta/\sqrt{x} tol ( x ) = 2 − 24 β / x S = Q × C S=Q\times C S = Q × C Q Q Q C = Toffolis + T / 4 C=\text{Toffolis}+T/4 C = Toffolis + T /4 T T T T T T n r o t n_{\mathrm{rot}} n rot 0.56 log 2 ( 1 / ϵ ) + 5.3 0.56\log_2(1/\epsilon)+5.3 0.56 log 2 ( 1/ ϵ ) + 5.3 [24 ] , where ϵ = π min x ( tol ( x ) − ∣ error ( x ) ∣ ) / n r o t \epsilon=\pi\min_x(\operatorname{tol}(x)-\lvert\operatorname{error}(x)\rvert)/n_{\mathrm{rot}} ϵ = π min x ( tol ( x ) − ∣ error ( x )∣) / n rot
Benchmark II: Reed-Solomon decoding for quantum advantage in optimal polynomial intersection
Problem introduction and importance
Decoded quantum interferometry (DQI) is a quantum algorithm that uses the quantum Fourier transform to reduce optimization problems to decoding problems [16 ] . We consider its application to the Optimal Polynomial Intersection (OPI) problem. An OPI instance specifies a set of allowed values F i ⊆ F q F_i \subseteq\mathbb{F}_q F i ⊆ F q m = q − 1 m=q-1 m = q − 1 α i \alpha_i α i i = 0 , … , m − 1 i=0,\ldots,m-1 i = 0 , … , m − 1 n < m n<m n < m Q x ( z ) = ∑ j = 0 n − 1 x j z j Q_x(z)=\sum_{j=0}^{n-1}x_jz^j Q x ( z ) = ∑ j = 0 n − 1 x j z j
f ( x ) = ∣ { i : Q x ( α i ) ∈ F i } ∣ . (32) f(x)=\left|\{i:Q_x(\alpha_i)\in F_i\}\right|. \tag*{(32)} f ( x ) = ∣ { i : Q x ( α i ) ∈ F i } ∣ . (32) Thus, a solution is a low-degree polynomial whose evaluations land in as many of the allowed sets as possible. The instances considered here have equally sized sets, ∣ F i ∣ = r |F_i|=r ∣ F i ∣ = r
OPI is a promising target for an early demonstration of classically verifiable quantum advantage. The quality of a sampled solution can be checked efficiently on a classical computer by evaluating f ( x ) f(x) f ( x ) [12 , 16 ] . OPI can be viewed as fitting a low-degree polynomial over a finite field to observations that specify several possible values at each input. Variants of this problem have been studied as Reed–Solomon list recovery in coding theory [25 ] and noisy polynomial reconstruction or interpolation in cryptography [26 , 27 ] .
Let B i j = α i j B_{ij}=\alpha_i^j B ij = α i j m × n m\times n m × n ( B x ) i = Q x ( α i ) (Bx)_i=Q_x(\alpha_i) ( B x ) i = Q x ( α i ) B x Bx B x C C C C ⊥ = ker ( B T ) C^\perp=\ker(B^T) C ⊥ = ker ( B T ) [ m , m − n , n + 1 ] [m,m-n,n+1] [ m , m − n , n + 1 ] d = n + 1 d=n+1 d = n + 1 ℓ \ell ℓ ∑ x P ( f ( x ) ) ∣ x ⟩ \sum_x P(f(x))|x\rangle ∑ x P ( f ( x )) ∣ x ⟩ P P P ℓ \ell ℓ Q x Q_x Q x ∣ P ( f ( x ) ) ∣ 2 |P(f(x))|^2 ∣ P ( f ( x )) ∣ 2 P P P [16 ] .
The connection between preparing ∑ x P ( f ( x ) ) ∣ x ⟩ \sum_x P(f(x))|x\rangle ∑ x P ( f ( x )) ∣ x ⟩ F q n \mathbb{F}_q^n F q n s , x ∈ F q n s,x\in\mathbb{F}_q^n s , x ∈ F q n s ⋅ x = ∑ j = 0 n − 1 s j x j ∈ F q s\cdot x=\sum_{j=0}^{n-1}s_jx_j\in\mathbb{F}_q s ⋅ x = ∑ j = 0 n − 1 s j x j ∈ F q
QFT ∣ s ⟩ = 1 q n ∑ x ∈ F q n χ ( s ⋅ x ) ∣ x ⟩ , (33) \operatorname{QFT}|s\rangle=\frac{1}{\sqrt{q^n}}\sum_{x\in\mathbb{F}_q^n}\chi(s\cdot x)|x\rangle, \tag*{(33)} QFT ∣ s ⟩ = q n 1 x ∈ F q n ∑ χ ( s ⋅ x ) ∣ x ⟩ , (33) where χ : F q → C \chi:\mathbb{F}_q\to\mathbb{C} χ : F q → C χ ( u + v ) = χ ( u ) χ ( v ) \chi(u+v)=\chi(u)\chi(v) χ ( u + v ) = χ ( u ) χ ( v ) q = 2 b q=2^b q = 2 b χ ( z ) = ( − 1 ) Tr ( z ) ∈ { ± 1 } \chi(z)=(-1)^{\operatorname{Tr}(z)}\in\{\pm1\} χ ( z ) = ( − 1 ) Tr ( z ) ∈ { ± 1 } QFT = QFT † \operatorname{QFT}=\operatorname{QFT}^{\dagger} QFT = QFT † h i ( z ) = 1 [ z ∈ F i ] h_i(z)=1[z\in F_i] h i ( z ) = 1 [ z ∈ F i ] b i ∈ F q n b_i\in\mathbb{F}_q^n b i ∈ F q n i i i B B B Q x ( α i ) = b i ⋅ x Q_x(\alpha_i)=b_i\cdot x Q x ( α i ) = b i ⋅ x f ( x ) = ∑ i = 0 m − 1 h i ( b i ⋅ x ) f(x)=\sum_{i=0}^{m-1}h_i(b_i\cdot x) f ( x ) = ∑ i = 0 m − 1 h i ( b i ⋅ x ) h i : F q → { 0 , 1 } h_i:\mathbb{F}_q\to\{0,1\} h i : F q → { 0 , 1 } (33) (n = 1 n=1 n = 1 h i ( z ) = ∑ y ∈ F q h i ^ ( y ) χ ( y z ) h_i(z)=\sum_{y\in\mathbb{F}_q}\widehat{h_i}(y)\chi(yz) h i ( z ) = ∑ y ∈ F q h i ( y ) χ ( y z ) z = b i ⋅ x z=b_i\cdot x z = b i ⋅ x y ∈ F q y\in\mathbb{F}_q y ∈ F q
h i ( b i ⋅ x ) = ∑ y ∈ F q h i ^ ( y ) χ ( y ( b i ⋅ x ) ) = ∑ y ∈ F q h i ^ ( y ) χ ( ( y b i ) ⋅ x ) . (34) h_i(b_i\cdot x)=\sum_{y\in\mathbb{F}_q}\widehat{h_i}(y)\chi\bigl(y(b_i\cdot x)\bigr)=\sum_{y\in\mathbb{F}_q}\widehat{h_i}(y)\chi\bigl((yb_i)\cdot x\bigr). \tag*{(34)} h i ( b i ⋅ x ) = y ∈ F q ∑ h i ( y ) χ ( y ( b i ⋅ x ) ) = y ∈ F q ∑ h i ( y ) χ ( ( y b i ) ⋅ x ) . (34) Because h i ( z ) ∈ { 0 , 1 } h_i(z)\in\{0,1\} h i ( z ) ∈ { 0 , 1 } h i k = h i h_i^k=h_i h i k = h i k ≥ 1 k\ge1 k ≥ 1 P ( f ( x ) ) = P ( ∑ i = 0 m − 1 h i ( b i ⋅ x ) ) P(f(x))=P\left(\sum_{i=0}^{m-1}h_i(b_i\cdot x)\right) P ( f ( x )) = P ( ∑ i = 0 m − 1 h i ( b i ⋅ x ) ) ∏ i ∈ S h i ( b i ⋅ x ) \prod_{i\in S}h_i(b_i\cdot x) ∏ i ∈ S h i ( b i ⋅ x ) S ⊆ { 0 , … , m − 1 } S\subseteq\{0,\ldots,m-1\} S ⊆ { 0 , … , m − 1 } ∣ S ∣ ≤ ℓ |S|\le\ell ∣ S ∣ ≤ ℓ (34) for each factor h i h_i h i y = e i ∈ F q y=e_i\in\mathbb{F}_q y = e i ∈ F q e ∈ F q m e\in\mathbb{F}_q^m e ∈ F q m e i = 0 e_i=0 e i = 0 i ∉ S i\notin S i ∈ / S e = ( e 0 , … , e m − 1 ) ∈ F q m e=(e_0,\ldots,e_{m-1})\in\mathbb{F}_q^m e = ( e 0 , … , e m − 1 ) ∈ F q m wt ( e ) ≤ ℓ \operatorname{wt}(e)\le\ell wt ( e ) ≤ ℓ χ ( u ) χ ( v ) = χ ( u + v ) \chi(u)\chi(v)=\chi(u+v) χ ( u ) χ ( v ) = χ ( u + v ) ∏ i ∈ S χ ( ( e i b i ) ⋅ x ) = χ ( ( ∑ i = 0 m − 1 e i b i ) ⋅ x ) = χ ( ( B T e ) ⋅ x ) \prod_{i\in S}\chi((e_ib_i)\cdot x)=\chi\left(\left(\sum_{i=0}^{m-1}e_ib_i\right)\cdot x\right)=\chi((B^Te)\cdot x) ∏ i ∈ S χ (( e i b i ) ⋅ x ) = χ ( ( ∑ i = 0 m − 1 e i b i ) ⋅ x ) = χ (( B T e ) ⋅ x ) (33) to each mode ∑ x χ ( ( B T e ) ⋅ x ) ∣ x ⟩ = q n QFT ∣ B T e ⟩ \sum_x\chi((B^Te)\cdot x)|x\rangle=\sqrt{q^n}\operatorname{QFT}|B^Te\rangle ∑ x χ (( B T e ) ⋅ x ) ∣ x ⟩ = q n QFT ∣ B T e ⟩
P ( f ( x ) ) = ∑ e ∈ F q m wt ( e ) ≤ ℓ a e χ ( B T e ⋅ x ) , so that ∑ x ∈ F q n P ( f ( x ) ) ∣ x ⟩ = q n QFT ( ∑ e ∈ F q m wt ( e ) ≤ ℓ a e ∣ B T e ⟩ ) . (35) P(f(x))=\sum_{\substack{e\in\mathbb{F}_q^m\\ \operatorname{wt}(e)\le\ell}}a_e\chi(B^Te\cdot x),\qquad\text{so that}\qquad\sum_{x\in\mathbb{F}_q^n}P(f(x))|x\rangle=\sqrt{q^n}\operatorname{QFT}\left(\sum_{\substack{e\in\mathbb{F}_q^m\\ \operatorname{wt}(e)\le\ell}}a_e|B^Te\rangle\right). \tag*{(35)} P ( f ( x )) = e ∈ F q m wt ( e ) ≤ ℓ ∑ a e χ ( B T e ⋅ x ) , so that x ∈ F q n ∑ P ( f ( x )) ∣ x ⟩ = q n QFT e ∈ F q m wt ( e ) ≤ ℓ ∑ a e ∣ B T e ⟩ . (35) This is why the polynomial degree ℓ \ell ℓ ∣ s ⟩ |s\rangle ∣ s ⟩ s = B T e s=B^Te s = B T e ℓ \ell ℓ
DQI prepares these components by constructing a superposition of sparse errors and computing their syndromes, giving ∑ e a e ∣ e ⟩ ∣ B T e ⟩ \sum_ea_e|e\rangle|B^Te\rangle ∑ e a e ∣ e ⟩ ∣ B T e ⟩ a e a_e a e P P P s = B T e s=B^Te s = B T e wt ( e ) ≤ ℓ \operatorname{wt}(e)\le\ell wt ( e ) ≤ ℓ C ⊥ C^\perp C ⊥ 2 ℓ ≤ n 2\ell\le n 2 ℓ ≤ n n n n n + 1 n+1 n + 1
∑ wt ( e ) ≤ ℓ a e ∣ e ⟩ ∣ B T e ⟩ ⟶ ∣ 0 ⟩ ∑ wt ( e ) ≤ ℓ a e ∣ B T e ⟩ . (36) \sum_{\substack{\operatorname{wt}(e)\le\ell}}a_e|e\rangle|B^Te\rangle\longrightarrow|0\rangle\sum_{\substack{\operatorname{wt}(e)\le\ell}}a_e|B^Te\rangle. \tag*{(36)} wt ( e ) ≤ ℓ ∑ a e ∣ e ⟩ ∣ B T e ⟩ ⟶ ∣0 ⟩ wt ( e ) ≤ ℓ ∑ a e ∣ B T e ⟩ . (36) An inverse Fourier transform on the syndrome register then produces the optimization state ∑ x P ( f ( x ) ∣ x ⟩ ∣ x ⟩ \sum_x P(f(x)|x\rangle|x\rangle ∑ x P ( f ( x ) ∣ x ⟩ ∣ x ⟩
The improved construction of Ref. [12 ] uses measurement-based uncomputation (MBUC) to generate and erase error symbols one at a time. For q = 2 b q=2^b q = 2 b b b b i i i e i e_i e i e i b i e_i b_i e i b i c i ∈ { 0 , 1 } b c_i\in\{0,1\}^b c i ∈ { 0 , 1 } b e i e_i e i ( − 1 ) c i ⋅ bin ( e i ) (-1)^{c_i\cdot\operatorname{bin}(e_i)} ( − 1 ) c i ⋅ bin ( e i ) m m m ( − 1 ) c ⋅ bin ( e ) (-1)^{c\cdot\operatorname{bin}(e)} ( − 1 ) c ⋅ bin ( e ) c = ( c 0 , … , c m − 1 ) c=(c_0,\ldots,c_{m-1}) c = ( c 0 , … , c m − 1 )
∣ ψ c ⟩ = ∑ wt ( e ) ≤ ℓ a e ( − 1 ) c ⋅ bin ( e ) ∣ B T e ⟩ ⟼ ∑ wt ( e ) ≤ ℓ a e ∣ B T e ⟩ . (37) |\psi_c\rangle=\sum_{\operatorname{wt}(e)\le\ell}a_e(-1)^{c\cdot\operatorname{bin}(e)}|B^Te\rangle\longmapsto\sum_{\operatorname{wt}(e)\le\ell}a_e|B^Te\rangle. \tag*{(37)} ∣ ψ c ⟩ = wt ( e ) ≤ ℓ ∑ a e ( − 1 ) c ⋅ bin ( e ) ∣ B T e ⟩ ⟼ wt ( e ) ≤ ℓ ∑ a e ∣ B T e ⟩ . (37) Here c ∈ { 0 , 1 } m b c\in\{0,1\}^{mb} c ∈ { 0 , 1 } mb bin ( e ) \operatorname{bin}(e) bin ( e ) s s s e ( s ) e(s) e ( s ) s s s c c c
The benchmark is to construct a kickmix circuit that implements this correction on computational-basis inputs:
U c ∣ s ⟩ ∣ 0 ⟩ W = ( − 1 ) c ⋅ bin ( e ( s ) ) ∣ s ⟩ ∣ 0 ⟩ W , s = B T e ( s ) , wt ( e ( s ) ) ≤ ℓ , (38) U_c|s\rangle|0\rangle_W=(-1)^{c\cdot\operatorname{bin}(e(s))}|s\rangle|0\rangle_W,\qquad s=B^Te(s),\quad\operatorname{wt}(e(s))\le\ell, \tag*{(38)} U c ∣ s ⟩ ∣0 ⟩ W = ( − 1 ) c ⋅ bin ( e ( s )) ∣ s ⟩ ∣0 ⟩ W , s = B T e ( s ) , wt ( e ( s )) ≤ ℓ , (38) Each test supplies a classical syndrome s s s ∣ s ⟩ |s\rangle ∣ s ⟩ c c c e ( s ) e(s) e ( s ) s s s W W W c c c U c U_c U c (37) by linearity. Implementations may recover error symbols one at a time or compute the required parity directly, allowing AI agents to optimize the decoding algorithm, arithmetic, and cleanup schedule within this contract.
Reference implementation
We consider an OPI instance over F 2 b \mathbb{F}_{2^b} F 2 b
( m , n , b , r ) = ( 4095 , 70 , 12 , 2016 ) , (39) (m,n,b,r)=(4095,70,12,2016), \tag*{(39)} ( m , n , b , r ) = ( 4095 , 70 , 12 , 2016 ) , (39) where r = ∣ F i ∣ = ( 2 b − 2 b / 2 ) / 2 r=|F_i|=(2^b-2^{b/2})/2 r = ∣ F i ∣ = ( 2 b − 2 b /2 ) /2 [ 4095 , 4025 , 71 ] [4095,4025,71] [ 4095 , 4025 , 71 ] F 2 12 \mathbb{F}_{2^{12}} F 2 12 ℓ = n / 2 = 35 \ell=n/2=35 ℓ = n /2 = 35 70 × 12 = 840 70\times12=840 70 × 12 = 840 c c c 4095 × 12 = 49 , 140 4095\times12=49,140 4095 × 12 = 49 , 140 F 2 12 ≅ F 2 [ α ] / ( α 12 + α 6 + α 4 + α + 1 ) \mathbb{F}_{2^{12}}\cong\mathbb{F}_2[\alpha]/(\alpha^{12}+\alpha^6+\alpha^4+\alpha+1) F 2 12 ≅ F 2 [ α ] / ( α 12 + α 6 + α 4 + α + 1 ) ( 1 , α , … , α 11 ) (1,\alpha,\ldots,\alpha^{11}) ( 1 , α , … , α 11 ) α j = α j \alpha_j=\alpha^j α j = α j j = 0 , … , m − 1 j=0,\ldots,m-1 j = 0 , … , m − 1 s i = ∑ j = 0 m − 1 e j α ( i + 1 ) j s_i=\sum_{j=0}^{m-1}e_j\alpha^{(i+1)j} s i = ∑ j = 0 m − 1 e j α ( i + 1 ) j i = 0 , … , n − 1 i=0,\ldots,n-1 i = 0 , … , n − 1
The baseline follows the decoder based on the Extended Euclidean Algorithm (EEA) in Ref. [12 ] . It has three main stages.
Solve the key equation. Interpret the syndrome as S ( z ) = ∑ j = 0 n − 1 s j z j S(z)=\sum_{j=0}^{n-1}s_jz^j S ( z ) = ∑ j = 0 n − 1 s j z j z 2 ℓ z^{2\ell} z 2 ℓ S ( z ) S(z) S ( z ) ℓ \ell ℓ σ \sigma σ Ω \Omega Ω
σ ( z ) S ( z ) ≡ Ω ( z ) ( m o d z 2 ℓ ) . (40) \sigma(z)S(z)\equiv\Omega(z)\pmod{z^{2\ell}}. \tag*{(40)} σ ( z ) S ( z ) ≡ Ω ( z ) ( mod z 2 ℓ ) . (40) The error-locator polynomial is proportional to ∏ i : e i ≠ 0 ( 1 − α i z ) \prod_{i:e_i\ne0}(1-\alpha_i z) ∏ i : e i = 0 ( 1 − α i z ) σ = 1 \sigma=1 σ = 1 Ω = 0 \Omega=0 Ω = 0 [12 , 28 ] , which shares registers between shrinking remainders and growing Bézout coefficients and stores temporary quotients in available register space. Its output gives explicit access to the Bézout coefficients needed by the next stage.
Recover symbols and correct their phases. A Chien search evaluates σ ( α i − 1 ) \sigma(\alpha_i^{-1}) σ ( α i − 1 ) m m m e i e_i e i Ω \Omega Ω σ ′ \sigma' σ ′ Z Z Z c i c_i c i ( − 1 ) c i ⋅ bin ( e i ) (-1)^{c_i\cdot\operatorname{bin}(e_i)} ( − 1 ) c i ⋅ bin ( e i )
Restore the syndrome and workspace. Uncompute the EEA from Step 1, restoring the original syndrome and returning all auxiliary registers to zero. This preserves the phase applied in Step 2, yielding Eq. (38) .
Our reference circuit stores polynomials in two arrays of n + 3 n+3 n + 3 b b b n n n O ( log n ) O(\log n) O ( log n ) 3 ( n + 1 ) + 9 = 222 3(n+1)+9=222 3 ( n + 1 ) + 9 = 222 F 2 12 \mathbb{F}_{2^{12}} F 2 12 α \alpha α F 2 \mathbb{F}_2 F 2 X X X C Z CZ C Z [12 ] .
In Step 2, each of the m = 4095 m=4095 m = 4095 σ \sigma σ σ ′ \sigma' σ ′ Ω \Omega Ω σ ′ ( α i − 1 ) \sigma'(\alpha_i^{-1}) σ ′ ( α i − 1 ) 2 ( n + 3 ) b + O ( log n ) 2(n+3)b+O(\log n) 2 ( n + 3 ) b + O ( log n ) 5.2 × 10 6 5.2\times10^6 5.2 × 1 0 6 [12 ] , which used a hand-tuned 51-Toffoli F 2 12 \mathbb{F}_{2^{12}} F 2 12 F 2 12 \mathbb{F}_{2^{12}} F 2 12
Verification strategy
DQI does not require a decoder that succeeds on every input. Every candidate polynomial produced by DQI is checked classically by evaluating f ( x ) f(x) f ( x ) (37) , each error e e e wt ( e ) ≤ ℓ \operatorname{wt}(e)\le\ell wt ( e ) ≤ ℓ ∣ a e ∣ 2 |a_e|^2 ∣ a e ∣ 2 s = B T e s=B^Te s = B T e ε \varepsilon ε s = B T e s=B^Te s = B T e e ∼ ∣ a e ∣ 2 e\sim|a_e|^2 e ∼ ∣ a e ∣ 2 c c c (3) and the fidelity argument after it, the output state of the decoder lies within trace distance 2 ε 2\sqrt{\varepsilon} 2 ε
Trace distance upper-bounds the total-variation distance between the outcome distributions of any measurement, so the implemented DQI circuit samples candidates within total-variation distance 2 ε 2\sqrt{\varepsilon} 2 ε p 0 ≥ 1 / 2 p_0\ge1/2 p 0 ≥ 1/2 ε \varepsilon ε
p ≥ 1 2 − 2 ε . (41) p\ge\frac{1}{2}-2\sqrt{\varepsilon}. \tag*{(41)} p ≥ 2 1 − 2 ε . (41) To test a submission, the verifier uses the seeded generator of Section II C to sample an error vector e ∼ ∣ a e ∣ 2 e\sim|a_e|^2 e ∼ ∣ a e ∣ 2 wt ( e ) ≤ ℓ \operatorname{wt}(e)\le\ell wt ( e ) ≤ ℓ ℓ \ell ℓ P P P [16 ] ), the symbol weight w = wt ( e ) ∈ { 0 , … , ℓ } w=\operatorname{wt}(e)\in\{0,\ldots,\ell\} w = wt ( e ) ∈ { 0 , … , ℓ } Pr [ wt ( e ) = w ] = v w \Pr[\operatorname{wt}(e)=w]=v_w Pr [ wt ( e ) = w ] = v w v = ( v 0 , … , v ℓ ) v=(v_0,\ldots,v_\ell) v = ( v 0 , … , v ℓ ) ( ℓ + 1 ) × ( ℓ + 1 ) (\ell+1)\times(\ell+1) ( ℓ + 1 ) × ( ℓ + 1 ) A A A
A w , w = w ( q − 2 r ) r ( q − r ) , A w − 1 , w = A w , w − 1 = w ( m − w + 1 ) . (42) A_{w,w}=\frac{w(q-2r)}{\sqrt{r(q-r)}},\qquad A_{w-1,w}=A_{w,w-1}=\sqrt{w(m-w+1)}. \tag*{(42)} A w , w = r ( q − r ) w ( q − 2 r ) , A w − 1 , w = A w , w − 1 = w ( m − w + 1 ) . (42) Conditioned on wt ( e ) = w \operatorname{wt}(e)=w wt ( e ) = w w w w { 0 , … , m − 1 } \{0,\ldots,m-1\} { 0 , … , m − 1 } e j e_j e j F q ∖ { 0 } \mathbb{F}_q\setminus\{0\} F q ∖ { 0 } E [ wt ( e ) ] ≈ 31.05 \mathbb{E}[\operatorname{wt}(e)]\approx31.05 E [ wt ( e )] ≈ 31.05 Pr [ wt ( e ) = 35 ] ≈ 2.8 % \Pr[\operatorname{wt}(e)=35]\approx2.8\% Pr [ wt ( e ) = 35 ] ≈ 2.8% n b nb nb s = B T e s=B^Te s = B T e m b mb mb c c c s s s c c c e e e ( − 1 ) e ⋅ b i n ( e ) (-1)^{e\cdot b_{\mathrm{in}}(e)} ( − 1 ) e ⋅ b in ( e ) k k k N = 102,400 N=102{,}400 N = 102 , 400 ε + \varepsilon^+ ε + ε \varepsilon ε 1 − 2 − 128 1-2^{-128} 1 − 2 − 128 (4) .
The score measures the cost to obtain an acceptable DQI candidate (lower is better). It weights qubits more heavily than Toffoli gates:
S = R ( ε + ) Q 2 T , R ( ε ) = ln ( 1 − p 0 ) ln ( 1 − ( p 0 − 2 ε ) ) = ln 2 − ln ( 1 2 + 2 ε ) . (43) S=R(\varepsilon^+)Q^2T,\qquad R(\varepsilon)=\frac{\ln(1-p_0)}{\ln(1-(p_0-2\sqrt{\varepsilon}))}=\frac{\ln2}{-\ln\left(\frac{1}{2}+2\sqrt{\varepsilon}\right)}. \tag*{(43)} S = R ( ε + ) Q 2 T , R ( ε ) = ln ( 1 − ( p 0 − 2 ε )) ln ( 1 − p 0 ) = − ln ( 2 1 + 2 ε ) ln 2 . (43) where we set p 0 = 1 / 2 p_0=1/2 p 0 = 1/2 Q Q Q T T T Q 2 T Q^2T Q 2 T
The factor R R R p p p t ( p ) = ln δ ′ / ln ( 1 − p ) t(p) = \ln\delta^{\prime}/\ln(1-p) t ( p ) = ln δ ′ / ln ( 1 − p ) δ ′ \delta^{\prime} δ ′ p 0 p_0 p 0 p 0 − 2 ε p_0-2\sqrt{\varepsilon} p 0 − 2 ε t ( p 0 − 2 ε ) / t ( p 0 ) = R ( ε ) t(p_0-2\sqrt{\varepsilon})/t(p_0)=R(\varepsilon) t ( p 0 − 2 ε ) / t ( p 0 ) = R ( ε ) δ ′ \delta^{\prime} δ ′ ε + = 10 − 5 \varepsilon^{+}=10^{-5} ε + = 1 0 − 5 10 − 4 10^{-4} 1 0 − 4 10 − 3 10^{-3} 1 0 − 3 10 − 2 10^{-2} 1 0 − 2 R ( ε + ) = 1.02 R(\varepsilon^{+})=1.02 R ( ε + ) = 1.02 1.06 1.06 1.06 1.21 1.21 1.21 1.94 1.94 1.94 ε + ≥ 1 / 16 \varepsilon^{+}\geq1/16 ε + ≥ 1/16 N = 102,400 N=102{,}400 N = 102 , 400 ε + ≈ 8.7 × 10 − 4 \varepsilon^{+}\approx8.7\times10^{-4} ε + ≈ 8.7 × 1 0 − 4 R ( ε + ) ≈ 1.19 R(\varepsilon^{+})\approx1.19 R ( ε + ) ≈ 1.19
Benchmark III: QROM and Majorana rotations for quantum simulations of chemistry
Problem introduction and importance
Quantum circuits for state-of-the-art ground-state energy estimation of a second-quantized electronic-structure Hamiltonian involve an important class of data-lookup circuits which are leading Toffoli cost. These costs are intimately related to the cost of specifying the Hamiltonian and are thus important primitives for algorithms encoding in second quantization. Thus their compilation is of interest beyond the benchmark provided here.
Key to many quantum algorithms is the specification of the operator being simulated. Our primary focus here will be using the linear combination of unitary approach which represents an operator as a sum of unitaries multiplied by a non-negative real weight. In the chemistry example considered the operator is the Hamiltonian and specifies how the system interacts. For a system of n n n O ( n 4 ) \mathcal{O}(n^{4}) O ( n 4 ) O ( n 4 ) \mathcal{O}(n^{4}) O ( n 4 )
Our goal will be to specify the Hamiltonian in Eq. (45) in the lowest space-time volume possible using an LCU where the unitaries are the rank-1 factorizations of the integral tensor. Almost all modern second quantized simulation techniques leverage compression of the two-body component of Eq. (45) which corresponds to writing the four tensor h p q r s ( 2 ) h_{pqrs}^{(2)} h pq r s ( 2 ) [11 ] and require the user to find optimal dimensions of the tensor network that achieve a good approximation of h p q r s ( 2 ) h_{pqrs}^{(2)} h pq r s ( 2 ) h p q ( 1 ) h_{pq}^{(1)} h pq ( 1 ) [11 ] as that requires a verification framework beyond kickmix, but we will use a proxy to ensure that users cannot artificially change the parameters to lower the space-time cost of the table lookups over the other components of the algorithm. The combination of classical optimization of the tensor network and how that leads to an efficient quantum circuit is representative of the type of large-scale computational workflow that informs the optimal circuit design.
Reference implementation
The block-encoding implementation is factored into self-contained components defined in this section. The user specifically optimizes over the following components that dominate the overall cost of the algorithm.
section III C 3, eq. (62) : The user explicitly optimizes the DFTHC factorization to Frobenius error ≤ ϵ D F T H C = 20 m H α \leq\epsilon_{\mathrm{DFTHC}}=20mH\alpha ≤ ϵ DFTHC = 20 m H α ( R , B , C ) (R,B,C) ( R , B , C ) u b ( ρ ) u_{b}^{(\rho)} u b ( ρ ) w b ( ρ ) w_{b}^{(\rho)} w b ( ρ ) ( S , β 1 ) (S,\beta_{1}) ( S , β 1 ) E g a p ( 2 Λ − E g a p ) \sqrt{E_{\mathrm{gap}}(2\Lambda-E_{\mathrm{gap}})} E gap ( 2Λ − E gap )
We provide the variational upper bound on ground state energy.
We provide the raw one- and two-electron integrals.
We provide a reference factorization that does not necessarily satisfy the accuracy requirements.
– We provide a reference address to bit-string map used to verify the output of the boolean oracles.
– We provide a reference tree.
The artifacts are input into the verification procedure of section III C 6 to verify a Frobenius error of at most ϵ D F T H C \epsilon_{\mathrm{DFTHC}} ϵ DFTHC
The following are provided to the user as black-boxes with heuristic cost functions. These are provided as black-boxes as there are too many possible implementations, and most of them would be difficult to verify. In particular, although many can be decomposed into Boolean oracle calls, these are highly sensitive to the particular implementation and are anyway subdominant in cost.
The overall score will be
S = λ e f f × T t o t × Q t o t , λ e f f = E g a p ( 2 Λ − E g a p ) , (44) S = \lambda_{\mathrm{eff}} \times T_{\mathrm{tot}} \times Q_{\mathrm{tot}}, \qquad\lambda_{\mathrm{eff}} = \sqrt{E_{\mathrm{gap}}(2\Lambda- E_{\mathrm{gap}})}, \tag*{(44)} S = λ eff × T tot × Q tot , λ eff = E gap ( 2Λ − E gap ) , (44) where T t o t T_{\mathrm{tot}} T tot k _k k N − 1 N-1 N − 1 Q t o t Q_{\mathrm{tot}} Q tot J b r o t J_{\mathrm{brot}} J brot k _k k k _k k < 1 < 1 < 1 ϵ D F T H C \epsilon_{\mathrm{DFTHC}} ϵ DFTHC k _k k Λ > 49.0 \Lambda> 49.0 Λ > 49.0 E g a p > 3.914898 E_{\mathrm{gap}} > 3.914898 E gap > 3.914898 λ e f f \lambda_{\mathrm{eff}} λ eff
Optimization of DFTHC factorization
Let [ N ] : = { 0 , … , N − 1 } [N] := \{0,\ldots,N-1\} [ N ] := { 0 , … , N − 1 }
H = ∑ p , q ∈ [ N ] h p q ( 1 ) E p q + 1 2 ∑ p , q , r , s ∈ [ N ] h p q r s ( 2 ) E p q E r s , E p q : = ∑ σ ∈ { 0 , 1 } a p σ † a q σ . (45) H = \sum_{p,q \in[N]} h^{(1)}_{pq} E_{pq} + \frac{1}{2}\sum_{p,q,r,s \in[N]} h^{(2)}_{pqrs} E_{pq}E_{rs}, \qquad E_{pq} := \sum_{\sigma\in\{0,1\}} a^\dagger_{p\sigma}a_{q\sigma}. \tag*{(45)} H = p , q ∈ [ N ] ∑ h pq ( 1 ) E pq + 2 1 p , q , r , s ∈ [ N ] ∑ h pq r s ( 2 ) E pq E r s , E pq := σ ∈ { 0 , 1 } ∑ a p σ † a q σ . (45) For our work we will search for representations of the Hamiltonian factorized into a form amenable to efficient circuit implementation. As such we can consider a related class of Hamiltonians that is spectrally equivalent in an η \eta η
H = ∑ p , q ∈ [ N ] ( h p q ( 1 ) − η 2 S p q + β 1 δ p q ) E p q + 1 2 ∑ p , q , r , s ∈ [ N ] ( h p q r s ( 2 ) + 1 2 ( S ⊗ I N + I N ⊗ S ) p q r s ) E p q E r s + ( − η β 1 ) I 4 N , (46) H = \sum_{p,q \in[N]}\left(h^{(1)}_{pq} - \frac{\eta}{2}S_{pq} + \beta_1\delta_{pq}\right)E_{pq} + \frac{1}{2}\sum_{p,q,r,s \in[N]}\left(h^{(2)}_{pqrs} + \frac{1}{2}(S \otimes I_N + I_N \otimes S)_{pqrs}\right)E_{pq}E_{rs} + (-\eta\beta_1)I_{4N}, \tag*{(46)} H = p , q ∈ [ N ] ∑ ( h pq ( 1 ) − 2 η S pq + β 1 δ pq ) E pq + 2 1 p , q , r , s ∈ [ N ] ∑ ( h pq r s ( 2 ) + 2 1 ( S ⊗ I N + I N ⊗ S ) pq r s ) E pq E r s + ( − η β 1 ) I 4 N , (46) for a symmetric matrix S S S β 1 \beta_1 β 1 ( S , β 1 ) (S,\beta_1) ( S , β 1 )
( h ( 1 ) , h ( 2 ) , R , B , C , H s y m , ϵ D F T H C ) , ϵ D F T H C : = 20 m H a = 2.0 × 10 − 2 H a , (47) \left(h^{(1)}, h^{(2)}, R, B, C, \mathcal{H}_{\mathrm{sym}}, \epsilon_{\mathrm{DFTHC}}\right), \qquad\epsilon_{\mathrm{DFTHC}} := 20\,\mathrm{mHa} = 2.0 \times10^{-2}\,\mathrm{Ha}, \tag*{(47)} ( h ( 1 ) , h ( 2 ) , R , B , C , H sym , ϵ DFTHC ) , ϵ DFTHC := 20 mHa = 2.0 × 1 0 − 2 Ha , (47) where H s y m \mathcal{H}_{\mathrm{sym}} H sym η \eta η R R R B B B C C C h p q r s ( 2 ) h^{(2)}_{pqrs} h pq r s ( 2 )
For every ρ ∈ [ R ] \rho\in[R] ρ ∈ [ R ] b ∈ [ B ] b \in[B] b ∈ [ B ]
u b ( ρ ) ∈ R N , ∥ u b ( ρ ) ∥ 2 = 1 , (48) u_b^{(\rho)} \in\mathbb{R}^N,\qquad\left\|u_b^{(\rho)}\right\|_2 = 1, \tag*{(48)} u b ( ρ ) ∈ R N , u b ( ρ ) 2 = 1 , (48) and, for every ( ρ , c ) ∈ [ R ] × [ C ] (\rho,c) \in[R] \times[C] ( ρ , c ) ∈ [ R ] × [ C ] w b ( ρ c ) w_b^{(\rho c)} w b ( ρ c ) W ( ρ c ) W^{(\rho c)} W ( ρ c )
X p q ( ρ c ) : = ∑ b ∈ [ B ] w b ( ρ c ) u b , p ( ρ ) u b , q ( ρ ) , h ~ p q r s ( 2 ) : = ∑ ρ ∈ [ R ] ∑ c ∈ [ C ] X p q ( ρ c ) X r s ( ρ c ) . (49) X_{pq}^{(\rho c)} := \sum_{b\in[B]} w_b^{(\rho c)}u_{b,p}^{(\rho)}u_{b,q}^{(\rho)},\qquad\widetilde{h}_{pqrs}^{(2)} := \sum_{\rho\in[R]}\sum_{c\in[C]} X_{pq}^{(\rho c)}X_{rs}^{(\rho c)}. \tag*{(49)} X pq ( ρ c ) := b ∈ [ B ] ∑ w b ( ρ c ) u b , p ( ρ ) u b , q ( ρ ) , h pq r s ( 2 ) := ρ ∈ [ R ] ∑ c ∈ [ C ] ∑ X pq ( ρ c ) X r s ( ρ c ) . (49) The corresponding DFTHC Hamiltonian is
H D F T H C ( Θ ) : = ∑ p , q h p q ( 1 ) ′ E p q + 1 2 ∑ ρ , c ( W ( ρ c ) I 4 N + ∑ p , q X p q ( ρ c ) E p q ) 2 − ( η β 1 + 1 2 ∑ ρ , c ∣ W ( ρ c ) ∣ 2 ) I 4 N , (50) H_{\mathrm{DFTHC}}(\Theta) := \sum_{p,q}h_{pq}^{(1)^\prime}E_{pq}+\frac{1}{2}\sum_{\rho,c}\left(W^{(\rho c)}I_{4^N}+\sum_{p,q}X_{pq}^{(\rho c)}E_{pq}\right)^2-\left(\eta\beta_1+\frac{1}{2}\sum_{\rho,c}\left|W^{(\rho c)}\right|^2\right)I_{4^N}, \tag*{(50)} H DFTHC ( Θ ) := p , q ∑ h pq ( 1 ) ′ E pq + 2 1 ρ , c ∑ ( W ( ρ c ) I 4 N + p , q ∑ X pq ( ρ c ) E pq ) 2 − ( η β 1 + 2 1 ρ , c ∑ W ( ρ c ) 2 ) I 4 N , (50) Θ : = { u b ( ρ ) , w b ( ρ c ) , W ( ρ c ) , h s y m ( 1 ) , β 1 } ρ , b , c (51) \Theta:= \left\{u_b^{(\rho)},w_b^{(\rho c)},W^{(\rho c)},h_{\mathrm{sym}}^{(1)},\beta_1\right\}_{\rho,b,c} \tag*{(51)} Θ := { u b ( ρ ) , w b ( ρ c ) , W ( ρ c ) , h sym ( 1 ) , β 1 } ρ , b , c (51) where h s y m ( 1 ) h_{\mathrm{sym}}^{(1)} h sym ( 1 ) β 1 \beta_1 β 1
h p q ( 1 ) ′ : = h p q ( 1 ) − η 2 S p q + β 1 δ p q − ∑ ρ , c W ( ρ c ) X p q ( ρ c ) . (52) h_{pq}^{(1)^\prime} := h_{pq}^{(1)}-\frac{\eta}{2}S_{pq}+\beta_1\delta_{pq}-\sum_{\rho,c}W^{(\rho c)}X_{pq}^{(\rho c)}. \tag*{(52)} h pq ( 1 ) ′ := h pq ( 1 ) − 2 η S pq + β 1 δ pq − ρ , c ∑ W ( ρ c ) X pq ( ρ c ) . (52) Consequently, expanding Eq. (50) gives
H D F T H C ( Θ ) = − η β 1 I 4 N + ∑ p , q ( h p q ( 1 ) − η 2 S p q + β 1 δ p q ) E p q + 1 2 ∑ p , q , r , s h ~ p q r s ( 2 ) E p q E r s . (53) H_{\mathrm{DFTHC}}(\Theta)=-\eta\beta_1I_{4^N}+\sum_{p,q}\left(h_{pq}^{(1)}-\frac{\eta}{2}S_{pq}+\beta_1\delta_{pq}\right)E_{pq}+\frac{1}{2}\sum_{p,q,r,s}\widetilde{h}_{pqrs}^{(2)}E_{pq}E_{rs}. \tag*{(53)} H DFTHC ( Θ ) = − η β 1 I 4 N + p , q ∑ ( h pq ( 1 ) − 2 η S pq + β 1 δ pq ) E pq + 2 1 p , q , r , s ∑ h pq r s ( 2 ) E pq E r s . (53) Thus the one-body interaction is represented exactly for every Θ \Theta Θ
ϵ f r o ( Θ ) : = ∥ h ( 2 ) − h ~ ( 2 ) ( Θ ) ∥ F = ( ∑ p , q , r , s ∣ h p q r s ( 2 ) − h ~ p q r s ( 2 ) ( Θ ) + 1 2 ( S ⊗ I N + I N ⊗ S ) p q r s ∣ 2 ) 1 / 2 . (54) \epsilon_{\mathrm{fro}}(\Theta):=\left\|h^{(2)}-\widetilde{h}^{(2)}(\Theta)\right\|_{\mathrm{F}}=\left(\sum_{p,q,r,s}\left|h_{pqrs}^{(2)}-\widetilde{h}_{pqrs}^{(2)}(\Theta)+\frac{1}{2}(S\otimes I_N+I_N\otimes S)_{pqrs}\right|^2\right)^{1/2}. \tag*{(54)} ϵ fro ( Θ ) := h ( 2 ) − h ( 2 ) ( Θ ) F = ( p , q , r , s ∑ h pq r s ( 2 ) − h pq r s ( 2 ) ( Θ ) + 2 1 ( S ⊗ I N + I N ⊗ S ) pq r s 2 ) 1/2 . (54) We will also check if the residual of the 1-electron component of the Hamiltonian is supplied correctly. To define the quantum-cost objective, let
w B ( ρ c ) : = W ( ρ c ) + ∑ b ∈ [ B ] w b ( ρ c ) . (55) w_B^{(\rho c)}:=W^{(\rho c)}+\sum_{b\in[B]}w_b^{(\rho c)}. \tag*{(55)} w B ( ρ c ) := W ( ρ c ) + b ∈ [ B ] ∑ w b ( ρ c ) . (55) The LCU one-norm, equivalently the normalization of the DFTHC Hamiltonian block encoding, is
Λ ( Θ ) : = ∥ h ( 1 ) ′ ∥ ∗ + 1 4 ∑ ρ , c ( ∣ w B ( ρ c ) ∣ + ∑ b ∈ [ B ] ∣ w b ( ρ c ) ∣ ) 2 , (56) \Lambda(\Theta):=\left\|h^{(1)^\prime}\right\|_*+\frac{1}{4}\sum_{\rho,c}\left(\left|w_B^{(\rho c)}\right|+\sum_{b\in[B]}\left|w_b^{(\rho c)}\right|\right)^2, \tag*{(56)} Λ ( Θ ) := h ( 1 ) ′ ∗ + 4 1 ρ , c ∑ w B ( ρ c ) + b ∈ [ B ] ∑ w b ( ρ c ) 2 , (56) where ∥ ⋅ ∥ ∗ \|\cdot\|_* ∥ ⋅ ∥ ∗
E S O S ( Θ ) : = − η β 1 − ( ∥ h ( 1 ) ′ ∥ ∗ − Tr h ( 1 ) ′ ) − 1 2 ∑ ρ , c ∣ W ( ρ c ) ∣ 2 , (57) E_{\mathrm{SOS}}(\Theta):=-\eta\beta_1-\left(\left\|h^{(1)^\prime}\right\|_*-\operatorname{Tr}h^{(1)^\prime}\right)-\frac{1}{2}\sum_{\rho,c}\left|W^{(\rho c)}\right|^2, \tag*{(57)} E SOS ( Θ ) := − η β 1 − ( h ( 1 ) ′ ∗ − Tr h ( 1 ) ′ ) − 2 1 ρ , c ∑ W ( ρ c ) 2 , (57) so that H D F T H C ( Θ ) − E S O S ( Θ ) I ⪰ 0 H_{\mathrm{DFTHC}}(\Theta)-E_{\mathrm{SOS}}(\Theta)I\succeq0 H DFTHC ( Θ ) − E SOS ( Θ ) I ⪰ 0
E g s ( Θ ) : = λ min ( H D F T H C ( Θ ) ∣ H s y m ) , (58) E_{\mathrm{gs}}(\Theta):=\lambda_{\min}\left(H_{\mathrm{DFTHC}}(\Theta)\big|_{\mathcal{H}_{\mathrm{sym}}}\right), \tag*{(58)} E gs ( Θ ) := λ m i n ( H DFTHC ( Θ ) H sym ) , (58) E g a p ( Θ ) : = E g s ( Θ ) − E S O S ( Θ ) , 0 ≤ E g a p ( Θ ) ≤ 2 Λ ( Θ ) . (59) E_{\mathrm{gap}}(\Theta):=E_{\mathrm{gs}}(\Theta)-E_{\mathrm{SOS}}(\Theta),\qquad0\leq E_{\mathrm{gap}}(\Theta)\leq2\Lambda(\Theta). \tag*{(59)} E gap ( Θ ) := E gs ( Θ ) − E SOS ( Θ ) , 0 ≤ E gap ( Θ ) ≤ 2Λ ( Θ ) . (59) Writing
δ g a p ( Θ ) : = E g a p ( Θ ) Λ ( Θ ) , g a m p ( Θ ) : = δ g a p ( Θ ) ( 2 − δ g a p ( Θ ) ) , (60) \delta_{\mathrm{gap}}(\Theta):=\frac{E_{\mathrm{gap}}(\Theta)}{\Lambda(\Theta)},\qquad g_{\mathrm{amp}}(\Theta):=\sqrt{\delta_{\mathrm{gap}}(\Theta)\left(2-\delta_{\mathrm{gap}}(\Theta)\right)}, \tag*{(60)} δ gap ( Θ ) := Λ ( Θ ) E gap ( Θ ) , g amp ( Θ ) := δ gap ( Θ ) ( 2 − δ gap ( Θ ) ) , (60) the spectrum-amplified block-encoding one-norm is the product
λ e f f ( Θ ) : = Λ ( Θ ) g a m p ( Θ ) = E g a p ( Θ ) ( 2 Λ ( Θ ) − E g a p ( Θ ) ) . (61) \lambda_{\mathrm{eff}}(\Theta) := \Lambda(\Theta)g_{\mathrm{amp}}(\Theta) = \sqrt{E_{\mathrm{gap}}(\Theta)\left(2\Lambda(\Theta)-E_{\mathrm{gap}}(\Theta)\right)}. \tag*{(61)} λ eff ( Θ ) := Λ ( Θ ) g amp ( Θ ) = E gap ( Θ ) ( 2Λ ( Θ ) − E gap ( Θ ) ) . (61) For fixed ( R , B , C ) (R,B,C) ( R , B , C )
Θ ⋆ ∈ argmin Θ λ e f f ( Θ ) subject to ϵ f r o ( Θ ) < ϵ D F T H C = 20 m H a , ∥ u b ( ρ ) ∥ 2 = 1 for all ρ ∈ [ R ] , b ∈ [ B ] . (62) \begin{aligned}
\Theta^\star&\in\underset{\Theta}{\operatorname{argmin}} && \lambda_{\mathrm{eff}}(\Theta)\\
&\text{subject to} && \epsilon_{\mathrm{fro}}(\Theta)<\epsilon_{\mathrm{DFTHC}}=20\ \mathrm{mHa},\\
&&& \lVert u_b^{(\rho)}\rVert_2=1\quad\text{for all }\rho\in[R],\ b\in[B]. \tag*{(62)}
\end{aligned} Θ ⋆ ∈ Θ argmin subject to λ eff ( Θ ) ϵ fro ( Θ ) < ϵ DFTHC = 20 mHa , ∥ u b ( ρ ) ∥ 2 = 1 for all ρ ∈ [ R ] , b ∈ [ B ] . (62) The output consists of Θ ⋆ \Theta^\star Θ ⋆ h ( 1 ) ′ ( Θ ⋆ ) h^{(1)\prime}(\Theta^\star) h ( 1 ) ′ ( Θ ⋆ ) ϵ f r o \epsilon_{\mathrm{fro}} ϵ fro Λ \Lambda Λ E S O S E_{\mathrm{SOS}} E SOS E g a p E_{\mathrm{gap}} E gap λ e f f \lambda_{\mathrm{eff}} λ eff ( R , B , C ) (R,B,C) ( R , B , C )
C B E ( R , B , C ) λ e f f ( Θ ⋆ ) , (63) \mathcal{C}_{\mathrm{BE}}(R,B,C)\lambda_{\mathrm{eff}}(\Theta^\star), \tag*{(63)} C BE ( R , B , C ) λ eff ( Θ ⋆ ) , (63) where C B E \mathcal{C}_{\mathrm{BE}} C BE
Verification of circuit for Boolean oracle
This is assumed to be a collection of K K K A N G L E S k \mathrm{ANGLES}_k ANGLES k ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ A N G L E S k \mathrm{ANGLES}_k ANGLES k ⌈ log 2 ( N + R B ) ⌉ \lceil\log_2(N+RB)\rceil ⌈ log 2 ( N + R B )⌉ a < N + R B a<N+RB a < N + R B N N N R B RB R B β \beta β β \beta β J J J [ 0 , 2 ) [0,2) [ 0 , 2 ) β − 1 \beta-1 β − 1 ∣ 0 ⟩ \lvert0\rangle ∣ 0 ⟩ π \pi π
Verification of Majorana rotations
Black-box routine: u ← TreeMajoranaCoefficients ( T , θ ⃗ ) \mathbf{u}\leftarrow\operatorname{TreeMajoranaCoefficients}(\mathcal{T},\vec{\theta}) u ← TreeMajoranaCoefficients ( T , θ )
Input: A rooted, child-ordered tree T \mathcal{T} T N N N θ ⃗ ∈ R N − 1 \vec{\theta}\in\mathbb{R}^{N-1} θ ∈ R N − 1 T \mathcal{T} T
Output: The coefficient vector u ( θ ⃗ ) ∈ R N \mathbf{u}(\vec{\theta})\in\mathbb{R}^{N} u ( θ ) ∈ R N x ∈ { 0 , 1 } x\in\{0,1\} x ∈ { 0 , 1 } γ u ( θ ⃗ ) , x : = ∑ v u v ( θ ⃗ ) γ v , x \gamma_{\mathbf{u}(\vec{\theta}),x}:=\sum_v u_v(\vec{\theta})\gamma_{v,x} γ u ( θ ) , x := ∑ v u v ( θ ) γ v , x U u ( θ ⃗ ) † γ 0 , x U u ( θ ⃗ ) = γ u ( θ ⃗ ) , x U_{\mathbf{u}(\vec{\theta})}^{\dagger}\gamma_{0,x}U_{\mathbf{u}(\vec{\theta})}=\gamma_{\mathbf{u}(\vec{\theta}),x} U u ( θ ) † γ 0 , x U u ( θ ) = γ u ( θ ) , x
Classical algorithm runtime: O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N )
Let T = ( V , E ) \mathcal{T}=(V,E) T = ( V , E ) V = { 0 , … , N − 1 } V=\{0,\ldots,N-1\} V = { 0 , … , N − 1 } 0 0 0 v v v γ v , x \gamma_{v,x} γ v , x x ∈ { 0 , 1 } x\in\{0,1\} x ∈ { 0 , 1 }
{ γ v , x , γ w , y } = 2 δ v w δ x y I . \{\gamma_{v,x},\gamma_{w,y}\}=2\delta_{vw}\delta_{xy}I. { γ v , x , γ w , y } = 2 δ v w δ x y I . Orient every edge away from the root, fix an order on the children of each vertex, and enumerate the edges in any parent-before-child order
e k = ( p k → q k ) , k = 1 , … , N − 1 , e_k=(p_k\to q_k),\qquad k=1,\ldots,N-1, e k = ( p k → q k ) , k = 1 , … , N − 1 , meaning that the edge entering a nonroot vertex precedes every edge leaving that vertex. Associate angle θ k \theta_k θ k e k e_k e k
For an oriented edge ( v → w ) (v\to w) ( v → w )
G v w ( θ ) : = exp [ θ 2 ( γ v , 0 γ w , 0 + γ v , 1 γ w , 1 ) ] . G_{vw}(\theta):=\exp\left[\frac{\theta}{2}(\gamma_{v,0}\gamma_{w,0}+\gamma_{v,1}\gamma_{w,1})\right]. G v w ( θ ) := exp [ 2 θ ( γ v , 0 γ w , 0 + γ v , 1 γ w , 1 ) ] . with convention
G v w ( θ ) † ( γ v , x γ w , x ) G v w ( θ ) = ( cos θ sin θ − sin θ cos θ ) ( γ v , x γ w , x ) , x ∈ { 0 , 1 } . G_{vw}(\theta)^\dagger
\begin{pmatrix}
\gamma_{v,x}\\
\gamma_{w,x}
\end{pmatrix}
G_{vw}(\theta)
=
\begin{pmatrix}
\cos\theta& \sin\theta\\
-\sin\theta& \cos\theta
\end{pmatrix}
\begin{pmatrix}
\gamma_{v,x}\\
\gamma_{w,x}
\end{pmatrix},
\qquad x\in\{0,1\}. G v w ( θ ) † ( γ v , x γ w , x ) G v w ( θ ) = ( cos θ − sin θ sin θ cos θ ) ( γ v , x γ w , x ) , x ∈ { 0 , 1 } . Set
U T ( θ ⃗ ) : = ∏ k = 1 N − 1 G p k q k ( θ k ) , U_T(\vec{\theta}):=\prod_{k=1}^{N-1}G_{p_kq_k}(\theta_k), U T ( θ ) := k = 1 ∏ N − 1 G p k q k ( θ k ) , with factors ordered from left to right in increasing k k k
Then there is a real unit vector u ( θ ⃗ ) ∈ R N u(\vec{\theta})\in\mathbb{R}^N u ( θ ) ∈ R N x = 0 , 1 x=0,1 x = 0 , 1
U T ( θ ⃗ ) † γ 0 , x U T ( θ ⃗ ) = ∑ v = 0 N − 1 u v ( θ ⃗ ) γ v , x = : γ u ( θ ⃗ ) , x , ∑ v = 0 N − 1 u v 2 = 1. (64) U_T(\vec{\theta})^\dagger\gamma_{0,x}U_T(\vec{\theta})
=\sum_{v=0}^{N-1}u_v(\vec{\theta})\gamma_{v,x}
=:\gamma_{u(\vec{\theta}),x},
\qquad
\sum_{v=0}^{N-1}u_v^2=1. \tag*{(64)} U T ( θ ) † γ 0 , x U T ( θ ) = v = 0 ∑ N − 1 u v ( θ ) γ v , x =: γ u ( θ ) , x , v = 0 ∑ N − 1 u v 2 = 1. (64) Thus a coherent implementation uses N − 1 N-1 N − 1
a. Classical coefficient algorithm. The coefficients in Eq. (64) are computed by applying the same planar rotations to e 0 \mathbf{e}_0 e 0
Input: ( p k , q k , θ k ) k = 1 N − 1 . Initialize: u ← ( 1 , 0 , … , 0 ) . For k = 1 , … , N − 1 : c ← cos θ k , s ← sin θ k , a ← u p k , b ← u q k , u p k ← c a − s b , u q k ← s a + c b . Return: u . (65) \begin{aligned}
\textbf{Input:}\quad& (p_k,q_k,\theta_k)_{k=1}^{N-1}.\\
\textbf{Initialize:}\quad& u\leftarrow(1,0,\ldots,0).\\
\textbf{For }k=1,\ldots,N-1:\quad& c\leftarrow\cos\theta_k,\quad s\leftarrow\sin\theta_k,\\
& a\leftarrow u_{p_k},\quad b\leftarrow u_{q_k},\\
& u_{p_k}\leftarrow ca-sb,\quad u_{q_k}\leftarrow sa+cb.\\
\textbf{Return:}\quad& u. \tag*{(65)}
\end{aligned} Input: Initialize: For k = 1 , … , N − 1 : Return: ( p k , q k , θ k ) k = 1 N − 1 . u ← ( 1 , 0 , … , 0 ) . c ← cos θ k , s ← sin θ k , a ← u p k , b ← u q k , u p k ← c a − s b , u q k ← s a + c b . u . (65) The ordering condition implies that q k q_k q k u q k = 0 u_{q_k}=0 u q k = 0 k k k
u q k ← u p k sin θ k , u p k ← u p k cos θ k . (66) u_{q_k}\leftarrow u_{p_k}\sin\theta_k,\qquad
u_{p_k}\leftarrow u_{p_k}\cos\theta_k. \tag*{(66)} u q k ← u p k sin θ k , u p k ← u p k cos θ k . (66) Each of the N − 1 N-1 N − 1 O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N ) (66) uses exactly 2 ( N − 1 ) 2(N-1) 2 ( N − 1 ) ∥ u ∥ 2 = 1 \lVert u\rVert_2=1 ∥ u ∥ 2 = 1
Equivalently, if the ordered children of v v v w 1 , … , w d v w_1,\ldots,w_{d_v} w 1 , … , w d v A v A_v A v v v v
A 0 = 1 , A w j = A v ( ∏ ℓ < j cos θ v w ℓ ) sin θ v w j , u v = A v ∏ ℓ = 1 d v cos θ v w ℓ . (67) A_0=1,\qquad
A_{w_j}=A_v\left(\prod_{\ell<j}\cos\theta_{vw_\ell}\right)\sin\theta_{vw_j},
\qquad
u_v=A_v\prod_{\ell=1}^{d_v}\cos\theta_{vw_\ell}. \tag*{(67)} A 0 = 1 , A w j = A v ℓ < j ∏ cos θ v w ℓ sin θ v w j , u v = A v ℓ = 1 ∏ d v cos θ v w ℓ . (67) For the path 0 → 1 → ⋯ → N − 1 0\to1\to\cdots\to N-1 0 → 1 → ⋯ → N − 1
Verification of DFTHC SOSSA block-encoding
Throughout, [ M ] : = { 0 , … , M − 1 } [M] := \{0,\ldots,M-1\} [ M ] := { 0 , … , M − 1 } 1 ≤ J ≤ N − 1 1 \le J \le N-1 1 ≤ J ≤ N − 1 K : = ⌈ ( N − 1 ) / J ⌉ K := \lceil(N-1)/J\rceil K := ⌈( N − 1 ) / J ⌉ T \mathcal{T} T Outer and Inner , and a family of Boolean angle oracles { A n g l e k } k = 0 K − 1 \{\mathrm{Angle}_k\}_{k=0}^{K-1} { Angle k } k = 0 K − 1
Let T = ( V , E ) \mathcal{T}=(V,E) T = ( V , E ) V = [ N ] V=[N] V = [ N ] e t = ( p t → q t ) e_t=(p_t\to q_t) e t = ( p t → q t ) t ∈ [ N − 1 ] t\in[N-1] t ∈ [ N − 1 ] a a a θ ⃗ a = ( θ a , 0 , … , θ a , N − 2 ) \vec{\theta}_a=(\theta_{a,0},\ldots,\theta_{a,N-2}) θ a = ( θ a , 0 , … , θ a , N − 2 ) { A n g l e k } k = 0 K − 1 \{\mathrm{Angle}_k\}_{k=0}^{K-1} { Angle k } k = 0 K − 1
u a : = TreeMajoranaCoefficients ( T , θ ⃗ a ) ∈ R N , U T , σ ( θ ⃗ a ) † γ 0 σ x U T , σ ( θ ⃗ a ) = ∑ p ∈ [ N ] u a , p γ p σ x . \mathbf{u}_a := \operatorname{TreeMajoranaCoefficients}(\mathcal{T},\vec{\theta}_a)\in\mathbb{R}^{N}, \qquad U_{\mathcal{T},\sigma}(\vec{\theta}_a)^\dagger\gamma_{0}\sigma_x U_{\mathcal{T},\sigma}(\vec{\theta}_a)=\sum_{p\in[N]}u_{a,p}\gamma_p\sigma_x. u a := TreeMajoranaCoefficients ( T , θ a ) ∈ R N , U T , σ ( θ a ) † γ 0 σ x U T , σ ( θ a ) = p ∈ [ N ] ∑ u a , p γ p σ x . The same real unit vector is obtained for both x ∈ { 0 , 1 } x\in\{0,1\} x ∈ { 0 , 1 } σ \sigma σ O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N ) (89) .
Let ℓ ∈ [ L ] \ell\in[L] ℓ ∈ [ L ] L = R C L=RC L = R C L L L R R R C C C B B B b ∈ [ B ] b\in[B] b ∈ [ B ] w ℓ b ∈ R w_{\ell b}\in\mathbb{R} w ℓ b ∈ R u ℓ b u_{\ell b} u ℓ b W ℓ W_\ell W ℓ w ℓ B = W ℓ + ∑ b ∈ [ B ] w ℓ b w_{\ell B}=W_\ell+\sum_{b\in[B]}w_{\ell b} w ℓ B = W ℓ + ∑ b ∈ [ B ] w ℓ b
X p q ( ℓ ) : = ∑ b ∈ [ B ] w ℓ b u ℓ b , p u ℓ b , q , h ~ p q r s ( 2 ) : = ∑ ℓ ∈ [ L ] X p q ( ℓ ) X r s ( ℓ ) , X_{pq}^{(\ell)}:=\sum_{b\in[B]}w_{\ell b}u_{\ell b,p}u_{\ell b,q}, \qquad\widetilde{h}_{pqrs}^{(2)}:=\sum_{\ell\in[L]}X_{pq}^{(\ell)}X_{rs}^{(\ell)}, X pq ( ℓ ) := b ∈ [ B ] ∑ w ℓ b u ℓ b , p u ℓ b , q , h pq r s ( 2 ) := ℓ ∈ [ L ] ∑ X pq ( ℓ ) X r s ( ℓ ) , h p q ( 1 ) ′ : = h p q ( 1 ) − η 2 S p q + β 1 δ p q − ∑ ℓ ∈ [ L ] W ℓ X p q ( ℓ ) . (68) h_{pq}^{(1)\prime}:=h_{pq}^{(1)}-\frac{\eta}{2}S_{pq}+\beta_1\delta_{pq}-\sum_{\ell\in[L]}W_\ell X_{pq}^{(\ell)}. \tag*{(68)} h pq ( 1 ) ′ := h pq ( 1 ) − 2 η S pq + β 1 δ pq − ℓ ∈ [ L ] ∑ W ℓ X pq ( ℓ ) . (68) Thus the one-electron tensor is absorbed exactly (though still verified); only h ~ ( 2 ) \widetilde{h}^{(2)} h ( 2 ) h ( 2 ) h^{(2)} h ( 2 )
Black-box input and block-encoding output
Inputs. Integers N N N R R R B B B C C C h ( 1 ) h^{(1)} h ( 1 ) h ( 2 ) h^{(2)} h ( 2 ) { w ℓ b , W ℓ } \{w_{\ell b},W_\ell\} { w ℓ b , W ℓ } T \mathcal{T} T b r o t b_{\mathrm{rot}} b rot 1 ≤ J ≤ N − 1 1\le J\le N-1 1 ≤ J ≤ N − 1
O u t e r ∣ 0 ⟩ = ∑ x ∈ X α x ∣ x ⟩ , Cost T ( O u t e r ) : = C o u t . (69) \mathrm{Outer}|0\rangle=\sum_{x\in\mathcal{X}}\alpha_x|x\rangle,\qquad\operatorname{Cost}_{\mathsf{T}}(\mathrm{Outer}):=C_{\mathrm{out}}. \tag*{(69)} Outer ∣0 ⟩ = x ∈ X ∑ α x ∣ x ⟩ , Cost T ( Outer ) := C out . (69) Here x x x
I n n e r : ∣ x ⟩ ∣ 0 ⟩ ⟼ ∣ x ⟩ ∑ y ∈ Y x β x , y ∣ y ⟩ , Cost T ( I n n e r ) : = C i n . (70) \mathrm{Inner}:|x\rangle|0\rangle\longmapsto|x\rangle\sum_{y\in\mathcal{Y}_x}\beta_{x,y}|y\rangle,\qquad\operatorname{Cost}_{\mathsf{T}}(\mathrm{Inner}):=C_{\mathrm{in}}. \tag*{(70)} Inner : ∣ x ⟩ ∣0 ⟩ ⟼ ∣ x ⟩ y ∈ Y x ∑ β x , y ∣ y ⟩ , Cost T ( Inner ) := C in . (70) Here y y y
For k ∈ [ K ] k\in[K] k ∈ [ K ]
ϑ a [ k ] : = ( ϑ a , k J , … , ϑ a , k J + J − 1 ) , \vartheta_a^{[k]}:=\left(\vartheta_{a,kJ},\ldots,\vartheta_{a,kJ+J-1}\right), ϑ a [ k ] := ( ϑ a , k J , … , ϑ a , k J + J − 1 ) , Δ ϑ a [ k ] : = { ϑ a [ 0 ] , k = 0 , ϑ a [ k ] ⊕ ϑ a [ k − 1 ] , k > 0. \Delta\vartheta_a^{[k]}:=
\begin{cases}
\vartheta_a^{[0]}, & k=0,\\
\vartheta_a^{[k]}\oplus\vartheta_a^{[k-1]}, & k>0.
\end{cases} Δ ϑ a [ k ] := { ϑ a [ 0 ] , ϑ a [ k ] ⊕ ϑ a [ k − 1 ] , k = 0 , k > 0. where ϑ a , t ∈ { 0 , 1 } b r o t \vartheta_{a,t}\in\{0,1\}^{b_{\mathrm{rot}}} ϑ a , t ∈ { 0 , 1 } b rot θ a , t / π \theta_{a,t}/\pi θ a , t / π t ≥ N − 1 t\ge N-1 t ≥ N − 1 k k k a a a
A n g l e k : ∣ a ⟩ ∣ z ⟩ ⟼ ∣ a ⟩ ∣ z ⊕ Δ ϑ a [ k ] ⟩ , Cost T ( A n g l e k ) : = C a n g , k ( J ) , k ∈ [ K ] . (71) \mathrm{Angle}_k:|a\rangle|z\rangle\longmapsto|a\rangle|z\oplus\Delta\vartheta_a^{[k]}\rangle,\qquad\operatorname{Cost}_{\mathsf{T}}(\mathrm{Angle}_k):=C_{\mathrm{ang},k}(J),\quad k\in[K]. \tag*{(71)} Angle k : ∣ a ⟩ ∣ z ⟩ ⟼ ∣ a ⟩ ∣ z ⊕ Δ ϑ a [ k ] ⟩ , Cost T ( Angle k ) := C ang , k ( J ) , k ∈ [ K ] . (71) XOR is componentwise over the J J J k > 0 k > 0 k > 0 ϑ a [ k − 1 ] \vartheta_a^{[k-1]} ϑ a [ k − 1 ] ϑ a [ k ] \vartheta_a^{[k]} ϑ a [ k ] C a n g p a s s ( J ) : = ∑ k ∈ [ K ] C a n g , k ( J ) C_{\mathrm{ang}}^{\mathrm{pass}}(J) := \sum_{k \in[K]} C_{\mathrm{ang},k}(J) C ang pass ( J ) := ∑ k ∈ [ K ] C ang , k ( J )
Output. A unitary U B E U_{\mathrm{BE}} U BE Π \Pi Π
Π U B E Π = ( H D F T H C − E S O S I Λ − I ) Π . (72) \Pi U_{\mathrm{BE}} \Pi= \left(\frac{H_{\mathrm{DFTHC}} - E_{\mathrm{SOS}} I}{\Lambda} - I\right)\Pi. \tag*{(72)} Π U BE Π = ( Λ H DFTHC − E SOS I − I ) Π. (72) The tree and these three oracle interfaces completely specify the implemented one- and two-electron block encoding. Each angle oracle is self-inverse, and each black-box inverse has the same assigned Toffoli cost.
Optimization scope. OUTER and INNER are never opened or re-synthesized: their contributions are represented only by heuristic values C o u t C_{\mathrm{out}} C out C i n C_{\mathrm{in}} C in T \mathcal{T} T { A n g l e k } k = 0 K − 1 \{\mathrm{Angle}_k\}_{k=0}^{K-1} { Angle k } k = 0 K − 1 C a n g , k ( J ) C_{\mathrm{ang},k}(J) C ang , k ( J ) C a n g p a s s ( J ) C_{\mathrm{ang}}^{\mathrm{pass}}(J) C ang pass ( J )
a. Streamed, batched, or parallel angles. Starting from z = 0 z = 0 z = 0 0 _0 0 k − 1 k - 1 k − 1 k _k k k k k 0 _0 0 z z z J = 1 J = 1 J = 1 J J J J J J J = N − 1 J = N - 1 J = N − 1 J J J
Verification and scoring
For FeMo-co (N = 76 N = 76 N = 76 η = 113 \eta= 113 η = 113 S = S 0 S = S_0 S = S 0 β 1 = β 1 , 0 \beta_1 = \beta_{1,0} β 1 = β 1 , 0
The verifier simulates every valid vector address, with 64 measurement samples per address, checking consistent angle words, unchanged addresses, clean workspace, and corrected phases during forward and reverse traversal. Each angle circuit is limited to 128 qubits. The recovered words determine the vectors used in all subsequent calculations.
Using these vectors, the verifier constructs X ( ρ c ) X(\rho^c) X ( ρ c ) h ( 1 ) ′ h^{(1)\prime} h ( 1 ) ′ h ( 1 ) ′ h^{(1)\prime} h ( 1 ) ′ h ~ ( 1 ) ′ \widetilde{h}^{(1)\prime} h ( 1 ) ′
ϵ 1 = ∥ h ~ ( 1 ) ′ − h ( 1 ) ′ ∥ F ≤ 0.020 H a , (73) \epsilon_1 = \left\lVert\widetilde{h}^{(1)\prime} - h^{(1)\prime} \right\rVert_{\mathrm{F}} \le0.020\ \mathrm{Ha}, \tag*{(73)} ϵ 1 = h ( 1 ) ′ − h ( 1 ) ′ F ≤ 0.020 Ha , (73) ϵ 2 = ∥ ∑ ρ , c X ( ρ c ) ⊗ X ( ρ c ) − h ( 2 ) − T ( S 0 ) ∥ F ≤ 0.020 H a , (74) \epsilon_2 = \left\lVert\sum_{\rho,c} X(\rho^c) \otimes X(\rho^c) - h^{(2)} - \mathcal{T}(S_0) \right\rVert_{\mathrm{F}} \le0.020\ \mathrm{Ha}, \tag*{(74)} ϵ 2 = ρ , c ∑ X ( ρ c ) ⊗ X ( ρ c ) − h ( 2 ) − T ( S 0 ) F ≤ 0.020 Ha , (74) where T ( S 0 ) = ( S 0 ⊗ I + I ⊗ S 0 ) / 2 \mathcal{T}(S_0) = (S_0 \otimes I + I \otimes S_0)/2 T ( S 0 ) = ( S 0 ⊗ I + I ⊗ S 0 ) /2
The verifier recomputes Λ \Lambda Λ E S O S E_{\mathrm{SOS}} E SOS
E g a p = E u p p e r − E S O S , (75) E_{\mathrm{gap}} = E_{\mathrm{upper}} - E_{\mathrm{SOS}}, \tag*{(75)} E gap = E upper − E SOS , (75) λ e f f = E g a p ( 2 Λ − E g a p ) , (76) \lambda_{\mathrm{eff}} = \sqrt{E_{\mathrm{gap}}(2\Lambda- E_{\mathrm{gap}})}, \tag*{(76)} λ eff = E gap ( 2Λ − E gap ) , (76) with the fixed, core-inclusive reference E u p p e r = − 22140.27335 H a E_{\mathrm{upper}} = -22140.27335\ \mathrm{Ha} E upper = − 22140.27335 Ha 0 < E g a p < 2 Λ 0 < E_{\mathrm{gap}} < 2\Lambda 0 < E gap < 2Λ E g a p > 3.914898 H a E_{\mathrm{gap}} > 3.914898\ \mathrm{Ha} E gap > 3.914898 Ha Λ > 49.0 H a \Lambda> 49.0\ \mathrm{Ha} Λ > 49.0 Ha
Accepted submissions receive the score
S = λ e f f T t o t Q t o t , (77) S = \lambda_{\mathrm{eff}} T_{\mathrm{tot}} Q_{\mathrm{tot}}, \tag*{(77)} S = λ eff T tot Q tot , (77) using simulated angle-circuit costs and the prescribed state-preparation and Givens-rotation cost model. Q t o t Q_{\mathrm{tot}} Q tot
The original baseline passes with unchanged score. These checks constrain coefficient errors and resource costs; they do not certify a 20 mHa many-electron energy error.
Classical quantities exposed to the user
Before evaluating either quantity below, the user decodes the angle words and applies TreeMajoranaCoefficients(T , θ ⃗ a \mathcal{T},\vec{\theta}_a T , θ a N + R B N+RB N + R B O ( ( N + R B ) N ) O((N+RB)N) O (( N + R B ) N ) O ( N ) O(N) O ( N )
a. Two-electron Frobenius error. The implemented error is
ϵ f r o : = ∥ h ( 2 ) − h ~ ( 2 ) ∥ F = [ ∑ p , q , r , s ( h p q r s ( 2 ) − ∑ ℓ X p q ( ℓ ) X r s ( ℓ ) ) 2 ] 1 / 2 . (78) \epsilon_{\mathrm{fro}} := \left\lVert h^{(2)}-\widetilde{h}^{(2)} \right\rVert_{\mathrm{F}} = \left[\sum_{p,q,r,s}\left(h^{(2)}_{pqrs}-\sum_{\ell}X^{(\ell)}_{pq}X^{(\ell)}_{rs}\right)^2\right]^{1/2}. \tag*{(78)} ϵ fro := h ( 2 ) − h ( 2 ) F = p , q , r , s ∑ ( h pq r s ( 2 ) − ℓ ∑ X pq ( ℓ ) X r s ( ℓ ) ) 2 1/2 . (78) It need not be evaluated by explicitly storing h ~ ( 2 ) \widetilde{h}^{(2)} h ( 2 ) H 2 ∈ R N 2 × N 2 \mathbf{H}_2\in\mathbb{R}^{N^2\times N^2} H 2 ∈ R N 2 × N 2 ( H 2 ) ( p q ) , ( r s ) = h p q r s ( 2 ) (\mathbf{H}_2)_{(pq),(rs)}=h^{(2)}_{pqrs} ( H 2 ) ( pq ) , ( r s ) = h pq r s ( 2 ) x ℓ = vec ( X ( ℓ ) ) x_\ell=\operatorname{vec}(X^{(\ell)}) x ℓ = vec ( X ( ℓ ) )
ϵ f r o 2 = ∥ h ( 2 ) ∥ F 2 − 2 ∑ ℓ x ℓ T H 2 x ℓ + ∑ ℓ , m ( x ℓ T x m ) 2 . (79) \epsilon_{\mathrm{fro}}^2=\left\lVert h^{(2)}\right\rVert_{\mathrm{F}}^2-2\sum_{\ell}x_\ell^{\mathsf{T}}\mathbf{H}_2x_\ell+\sum_{\ell,m}\left(x_\ell^{\mathsf{T}}x_m\right)^2. \tag*{(79)} ϵ fro 2 = h ( 2 ) F 2 − 2 ℓ ∑ x ℓ T H 2 x ℓ + ℓ , m ∑ ( x ℓ T x m ) 2 . (79) With ordinary dense floating-point contractions this costs
O ( N 4 + L B N 2 + L N 4 + L 2 N 2 ) (80) O(N^4+LBN^2+LN^4+L^2N^2) \tag*{(80)} O ( N 4 + L B N 2 + L N 4 + L 2 N 2 ) (80) operations, with dominant contraction cost O ( L N 4 ) = O ( R C N 4 ) O(LN^4)=O(RCN^4) O ( L N 4 ) = O ( R C N 4 ) B , C = Θ ( N ) B,C=\Theta(N) B , C = Θ ( N ) R = O ~ ( 1 ) R=\widetilde{O}(1) R = O ( 1 ) ϵ f r o \epsilon_{\mathrm{fro}} ϵ fro O ~ ( N 5 ) \widetilde{O}(N^5) O ( N 5 ) O ( N 5 ) O(N^5) O ( N 5 ) R R R
b. 1-electron Frobenius norm error We also check the loss with respect to the correct 1-electron coefficients. This is computed as a function of the tensors supplied by the user S S S W ρ , c W^{\rho,c} W ρ , c X p q ρ c X^{\rho c}_{pq} X pq ρ c H D F T H C H_{\mathrm{DFTHC}} H DFTHC N N N
h p q ( 1 ) ′ = h p q ( 1 ) − η 2 S p q + β 1 δ p q − ∑ ℓ ∈ [ L ] W ℓ X p q ( ℓ ) (81) h^{(1)\prime}_{pq}=h^{(1)}_{pq}-\frac{\eta}{2}S_{pq}+\beta_1\delta_{pq}-\sum_{\ell\in[L]}W_\ell X^{(\ell)}_{pq} \tag*{(81)} h pq ( 1 ) ′ = h pq ( 1 ) − 2 η S pq + β 1 δ pq − ℓ ∈ [ L ] ∑ W ℓ X pq ( ℓ ) (81) where X p q ( ℓ ) X^{(\ell)}_{pq} X pq ( ℓ ) R B RB R B
ϵ 1 , f r o = ∥ h p q ( 1 ) − h ~ p q ( 1 ) ∥ (82) \epsilon_{1,\mathrm{fro}}=\left\lVert h^{(1)}_{pq}-\widetilde{h}^{(1)}_{pq}\right\rVert \tag*{(82)} ϵ 1 , fro = h pq ( 1 ) − h pq ( 1 ) (82) is small and below the given threshold of 20 mHa. To do this we first construct h ( 1 ) ′ h^{(1)\prime} h ( 1 ) ′
h ~ ( 1 ) ′ = ∑ j λ j v ~ j v ~ j † (83) \widetilde{h}^{(1)\prime}=\sum_j\lambda_j\widetilde{v}_j\widetilde{v}_j^\dagger \tag*{(83)} h ( 1 ) ′ = j ∑ λ j v j v j † (83) where λ j \lambda_j λ j v ~ j \widetilde{v}_j v j β \beta β h ~ ( 1 ) \widetilde{h}^{(1)} h ( 1 )
h ~ p q ( 1 ) = h ~ p q ( 1 ) ′ + η 2 S p q − β 1 δ p q + ∑ ρ , c W ρ , c X ~ p q ρ c (84) \widetilde{h}^{(1)}_{pq}=\widetilde{h}^{(1)\prime}_{pq}+\frac{\eta}{2}S_{pq}-\beta_1\delta_{pq}+\sum_{\rho,c}W^{\rho,c}\widetilde{X}^{\rho c}_{pq} \tag*{(84)} h pq ( 1 ) = h pq ( 1 ) ′ + 2 η S pq − β 1 δ pq + ρ , c ∑ W ρ , c X pq ρ c (84) which enables the ϵ 1 , f r o \epsilon_{1,\mathrm{fro}} ϵ 1 , fro
c. Block-encoding normalization. Diagonalize the exactly corrected one-body matrix as
h ( 1 ) ′ = ∑ j w + , j u + , j u + , j T − ∑ j w − , j u − , j u − , j T , w ± , j ≥ 0. (85) h^{(1)\prime}=\sum_j w_{+,j}u_{+,j}u_{+,j}^{\mathsf{T}}-\sum_j w_{-,j}u_{-,j}u_{-,j}^{\mathsf{T}},\qquad w_{\pm,j}\geq0. \tag*{(85)} h ( 1 ) ′ = j ∑ w + , j u + , j u + , j T − j ∑ w − , j u − , j u − , j T , w ± , j ≥ 0. (85) The normalization is the sum of its one- and two-electron contributions:
Λ = ∑ j ( w + , j + w − , j ) ⏟ Λ 1 = ∥ h ( 1 ) ∥ ∗ + 1 4 ∑ ℓ ∈ [ L ] ( ∣ w ℓ B ∣ + ∑ b ∈ [ B ] ∣ w ℓ b ∣ ) 2 ⏟ Λ 2 . (86) \Lambda= \underbrace{\sum_{j} (w_{+,j} + w_{-,j})}_{\Lambda_{1} = \lVert h^{(1)}\rVert_*} + \underbrace{\frac{1}{4}\sum_{\ell\in[L]} \left( \lvert w_{\ell B}\rvert+ \sum_{b \in[B]} \lvert w_{\ell b}\rvert\right)^2}_{\Lambda_{2}}. \tag*{(86)} Λ = Λ 1 = ∥ h ( 1 ) ∥ ∗ j ∑ ( w + , j + w − , j ) + Λ 2 4 1 ℓ ∈ [ L ] ∑ ∣ w ℓ B ∣ + b ∈ [ B ] ∑ ∣ w ℓ b ∣ 2 . (86) Equations (78) and (86) are the only Hamiltonian-quality quantities required by the user-facing optimizer.
Outlook
The benchmarks in this paper evaluate flat gate streams using logical qubit and Toffoli counts. While this setup gives AI agents a fast inner loop for subroutine optimization, moving toward end-to-end fault-tolerant compilation suggests two main directions:
Compiling directly to a target QEC architecture: Logical qubit and Toffoli counts are coarse proxies for physical hardware cost. In practice, hardware capabilities dictate the choice of quantum error-correcting (QEC) code. That code choice then determines how the architecture implements logical gates and what overheads arise from magic-state distillation, qubit routing, and reaction time [7 , 10 , 12 , 29 , 30 ] . Future benchmarks can keep kickmix simulation to check logical correctness while scoring submissions with a QEC cost model for a target architecture. AI agents can then co-optimize the logical arithmetic, qubit layout, and gate schedule to minimize physical qubits and runtime.
Machine-checkable proofs of circuit correctness: Classical simulation applies only to kickmix circuits, and a flat .kmb gate stream only tests a single problem size. Many quantum subroutines, such as state preparation, quantum signal processing, and fermionic rotations, create superpositions internally and cannot be checked by a kickmix simulator. When human researchers optimize these subroutines, they design parameterized constructions for any input size n n n [31 ] mature, submissions could pair a circuit generator with a proof of its correctness and cost. Because a proof assistant checks algebraic identities rather than simulating basis states, it can certify arbitrary quantum unitaries regardless of Hilbert space dimension. This shift would let AI agents optimize general quantum algorithms end to end and discover parameterized constructions that hold across all problem sizes.
Together, fast verifiers, architecture-aware cost models, and formal proofs can democratize fault-tolerant compilation, opening quantum algorithm design to researchers and engineers across all fields.
Acknowledgments
While we were finalizing this manuscript and preparing the competition, Ref. [32 ] appeared on the arXiv. That work also builds on the kickmix model of Ref. [13 ] to optimize fault-tolerant quantum subroutines with AI agents, and studies many of the same benchmark problems, including Reed–Solomon syndrome decoding for DQI, Coulomb-potential evaluation, and quantum lookup tables for electronic structure. Although there is significant overlap between the two studies, this work was developed independently and in parallel. We note that we also spoke about accelerating these benchmarks two weeks prior in a recorded talk at KTH [33 ] . Our internal results show that all three benchmarks can be further improved beyond the circuits reported to date, so the competition remains open and interesting. Taken together, the strong results in Ref. [32 ] underscore the effectiveness of pairing fast classical verification with AI-driven search, highlighting both the importance of this approach and the relevance of the competition.
References
low-overhead transversal architectures for reconfigurable atom arrays (2025).
[30 ] M. Cain, Q. Xu, R. King, L. R. B. Picard, H. Levine, M. Endres, J. Preskill, H.-Y. Huang, and D. Bluvstein, Shor’s algorithm is possible with as few as 10,000 reconfigurable atomic qubits (2026).
[31 ] L. de Moura, S. Kong, J. Avigad, F. van Doorn, and J. von Raumer, The Lean theorem prover (system description), in Automated Deduction - CADE-25 (Springer International Publishing, 2015) pp. 378–388.
[32 ] A. Gu, A. H. Kemeny, S. F. Yelin, and M. D. Lukin, Discovering state-of-the-art quantum circuits with large language models (2026), arXiv:2610.22367 [quant-ph].
[33 ] R. Babbush, Accelerating the field of quantum algorithms with AI, Recorded talk at the KITP Conference: Artificial Intelligence at the Quantum Frontier (2026).
DFTHC block encoding circuits
This appendix expands the tree rotation and the opaque addresses in section III C 6 into the labels and select logic used by the DFTHC block encoding. None of these labels is an optimization variable in the user-facing problem.
Tree Majorana coefficients and basis rotation
Write the ordered edges of T \mathcal{T} T e t = ( p t → q t ) e_t=(p_t\to q_t) e t = ( p t → q t ) t ∈ [ N − 1 ] t\in[N-1] t ∈ [ N − 1 ] σ \sigma σ
G p q ( σ ) ( θ ) : = exp [ θ 2 ( γ p σ 0 γ q σ 0 + γ p σ 1 γ q σ 1 ) ] , (87) G_{p q}^{(\sigma)}(\theta):=\exp\left[\frac{\theta}{2}\left(\gamma_{p\sigma0}\gamma_{q\sigma0}+\gamma_{p\sigma1}\gamma_{q\sigma1}\right)\right], \tag*{(87)} G pq ( σ ) ( θ ) := exp [ 2 θ ( γ p σ 0 γ q σ 0 + γ p σ 1 γ q σ 1 ) ] , (87) and let
U T , σ ( θ ⃗ ) : = ∏ t = 0 N − 2 G p t q t ( σ ) ( θ t ) , (88) U_{\mathcal{T},\sigma}(\vec{\theta}):=\prod_{t=0}^{N-2}G_{p_tq_t}^{(\sigma)}(\theta_t), \tag*{(88)} U T , σ ( θ ) := t = 0 ∏ N − 2 G p t q t ( σ ) ( θ t ) , (88) with factors written from left to right in increasing t t t
u ← ( 1 , 0 , … , 0 ) , \mathbf{u}\leftarrow(1,0,\ldots,0), u ← ( 1 , 0 , … , 0 ) , ( u p t , u q t ) ← ( cos θ t u p t − sin θ t u q t , sin θ t u p t + cos θ t u q t ) , t = 0 , … , N − 2. (89) (u_{p_t},u_{q_t})\leftarrow(\cos\theta_t\,u_{p_t}-\sin\theta_t\,u_{q_t},\sin\theta_t\,u_{p_t}+\cos\theta_t\,u_{q_t}),\qquad t=0,\ldots,N-2. \tag*{(89)} ( u p t , u q t ) ← ( cos θ t u p t − sin θ t u q t , sin θ t u p t + cos θ t u q t ) , t = 0 , … , N − 2. (89) It returns a unit vector and satisfies, simultaneously for x = 0 , 1 x=0,1 x = 0 , 1
U T , σ ( θ ⃗ ) † γ 0 σ x U T , σ ( θ ⃗ ) = ∑ p ∈ [ N ] u p γ p σ x . (90) U_{\mathcal{T},\sigma}(\vec{\theta})^\dagger\gamma_{0\sigma x}U_{\mathcal{T},\sigma}(\vec{\theta})=\sum_{p\in[N]}u_p\gamma_{p\sigma x}. \tag*{(90)} U T , σ ( θ ) † γ 0 σ x U T , σ ( θ ) = p ∈ [ N ] ∑ u p γ p σ x . (90) Each of the N − 1 N-1 N − 1 O ( N ) O(N) O ( N ) O ( N ) O(N) O ( N ) 0 _0 0 k _k k
The quantum circuit implementing the t t h t^{\mathrm{th}} t th b t b_t b t 2 b t − 1 2b_t-1 2 b t − 1
Sum-of-squares generators
Let r ∈ [ R ] r\in[R] r ∈ [ R ] c ∈ [ C ] c\in[C] c ∈ [ C ] b ∈ [ B ] b\in[B] b ∈ [ B ] σ ∈ { 0 , 1 } \sigma\in\{0,1\} σ ∈ { 0 , 1 } G ∈ { D 1 , Q 1 , S F } G\in\{\mathrm{D}_1,\mathrm{Q}_1,\mathrm{SF}\} G ∈ { D 1 , Q 1 , SF } u u u
γ u σ x : = ∑ p = 0 N − 1 u p γ p σ x , x ∈ { 0 , 1 } . (91) \gamma_{u\sigma x}:=\sum_{p=0}^{N-1}u_p\gamma_{p\sigma x},\qquad x\in\{0,1\}. \tag*{(91)} γ u σ x := p = 0 ∑ N − 1 u p γ p σ x , x ∈ { 0 , 1 } . (91) The one-electron positive and negative eigenspaces and the two-electron squares are represented by
O D 1 , j : = w + , j γ u + , j σ 0 + i γ u + , j σ 1 2 , (92) O_{\mathrm{D}_1,j}:=\sqrt{w_{+,j}}\frac{\gamma_{u_{+,j}}\sigma_0+i\gamma_{u_{+,j}}\sigma_1}{2}, \tag*{(92)} O D 1 , j := w + , j 2 γ u + , j σ 0 + i γ u + , j σ 1 , (92) O Q 1 , j : = w − , j γ u − , j σ 0 − i γ u − , j σ 1 2 , (93) O_{\mathrm{Q}_1,j}:=\sqrt{w_{-,j}}\frac{\gamma_{u_{-,j}}\sigma_0-i\gamma_{u_{-,j}}\sigma_1}{2}, \tag*{(93)} O Q 1 , j := w − , j 2 γ u − , j σ 0 − i γ u − , j σ 1 , (93) O S F , r c : = 1 2 ( w B ( r c ) I + i 2 ∑ b ∈ [ B ] ∑ σ ∈ { 0 , 1 } w b ( r c ) γ u b ( r ) σ 0 γ u b ( r ) σ 1 ) , (94) O_{\mathrm{SF},rc}:=\frac{1}{\sqrt{2}}\left(w_B^{(rc)}I+\frac{i}{2}\sum_{b\in[B]}\sum_{\sigma\in\{0,1\}}w_b^{(rc)}\gamma_{u_b^{(r)}}\sigma_0\gamma_{u_b^{(r)}}\sigma_1\right), \tag*{(94)} O SF , r c := 2 1 w B ( r c ) I + 2 i b ∈ [ B ] ∑ σ ∈ { 0 , 1 } ∑ w b ( r c ) γ u b ( r ) σ 0 γ u b ( r ) σ 1 , (94) where w B ( r c ) = W ( r c ) + ∑ b w b ( r c ) w_B^{(rc)}=W^{(rc)}+\sum_b w_b^{(rc)} w B ( r c ) = W ( r c ) + ∑ b w b ( r c )
H D F T H C − E S O S I = ∑ j , σ ( O D 1 , j † O D 1 , j + O Q 1 , j † O Q 1 , j ) + ∑ r , c O S F , r c † O S F , r c , (95) H_{\mathrm{DFTHC}}-E_{\mathrm{SOS}}I=\sum_{j,\sigma}\left(O_{\mathrm{D}_1,j}^{\dagger}O_{\mathrm{D}_1,j}+O_{\mathrm{Q}_1,j}^{\dagger}O_{\mathrm{Q}_1,j}\right)+\sum_{r,c}O_{\mathrm{SF},rc}^{\dagger}O_{\mathrm{SF},rc}, \tag*{(95)} H DFTHC − E SOS I = j , σ ∑ ( O D 1 , j † O D 1 , j + O Q 1 , j † O Q 1 , j ) + r , c ∑ O SF , r c † O SF , r c , (95) with
E S O S = − 2 ∑ j w − , j − 1 2 ∑ r , c ∣ W ( r c ) ∣ 2 . (96) E_{\mathrm{SOS}}=-2\sum_j w_{-,j}-\frac{1}{2}\sum_{r,c}\left|W^{(rc)}\right|^2. \tag*{(96)} E SOS = − 2 j ∑ w − , j − 2 1 r , c ∑ W ( r c ) 2 . (96) The generator normalizations are
λ D 1 , j = w + , j , λ Q 1 , j = w − , j , λ S F , r c = 1 2 ( ∣ w B ( r c ) ∣ + ∑ b ∣ w b ( r c ) ∣ ) , (97) \lambda_{\mathrm{D}_1,j}=\sqrt{w_{+,j}},\qquad\lambda_{\mathrm{Q}_1,j}=\sqrt{w_{-,j}},\qquad\lambda_{\mathrm{SF},rc}=\frac{1}{\sqrt{2}}\left(\left|w_B^{(rc)}\right|+\sum_b\left|w_b^{(rc)}\right|\right), \tag*{(97)} λ D 1 , j = w + , j , λ Q 1 , j = w − , j , λ SF , r c = 2 1 ( w B ( r c ) + b ∑ w b ( r c ) ) , (97) and Λ \Lambda Λ (86) equals one half of the sum of their squared normalizations, including both spins.
Registers, address maps, and select bits
Let N + N_+ N + N − N_- N − h ( 1 ) ′ h^{(1)'} h ( 1 ) ′ N + + N − = N N_++N_-=N N + + N − = N x o ∈ [ N + R C ] x_o\in[N+RC] x o ∈ [ N + R C ]
x o ( G , r , c ) = { r , G = D 1 , r ∈ [ N + ] , N + + r , G = Q 1 , r ∈ [ N − ] , N + r C + c , G = S F , r ∈ [ R ] , c ∈ [ C ] . (98) x_o(G,r,c)=
\begin{cases}
r, & G=\mathrm{D}_1,\quad r\in[N_+],\\
N_++r, & G=\mathrm{Q}_1,\quad r\in[N_-],\\
N+rC+c, & G=\mathrm{SF},\quad r\in[R],\ c\in[C].
\end{cases} \tag*{(98)} x o ( G , r , c ) = ⎩ ⎨ ⎧ r , N + + r , N + r C + c , G = D 1 , r ∈ [ N + ] , G = Q 1 , r ∈ [ N − ] , G = SF , r ∈ [ R ] , c ∈ [ C ] . (98) For a two-electron generator, b < B b<B b < B b = B b=B b = B μ ∈ [ N + R B ] \mu\in[N+RB] μ ∈ [ N + R B ] N N N R B RB R B u b ( r ) u_b^{(r)} u b ( r ) a a a (71) ; the batch label k k k A N G L E k \mathrm{ANGLE}_k ANGLE k
The remaining finite select data comprise: the generator type G G G x ∈ { 0 , 1 } x\in\{0,1\} x ∈ { 0 , 1 } b = B b=B b = B
Circuit order and cost accounting
At a high level, one generator is block-encoded by
OUTER ⟶ [ decode ⟶ INNER ⟶ spin select ] ⟶ { ANGLE k ⟶ Givens batch k } k = 0 K − 1 \operatorname{OUTER}\longrightarrow[\operatorname{decode}\longrightarrow\operatorname{INNER}\longrightarrow\operatorname{spin\ select}]\longrightarrow\{\operatorname{ANGLE}_k\longrightarrow\operatorname{Givens\ batch\ }k\}_{k=0}^{K-1} OUTER ⟶ [ decode ⟶ INNER ⟶ spin select ] ⟶ { ANGLE k ⟶ Givens batch k } k = 0 K − 1 ⟶ Majorana select ⟶ { inverse Givens batch k ⟶ ANGLE k } k = K − 1 0 \longrightarrow\operatorname{Majorana\ select}\longrightarrow\{\operatorname{inverse\ Givens\ batch\ }k\longrightarrow\operatorname{ANGLE}_k\}_{k=K-1}^{0} ⟶ Majorana select ⟶ { inverse Givens batch k ⟶ ANGLE k } k = K − 1 0 ⟶ [ INNER † ⟶ reflect and square ] ⟶ OUTER † . (99) \longrightarrow[\operatorname{INNER}^{\dagger}\longrightarrow\operatorname{reflect\ and\ square}]\longrightarrow\operatorname{OUTER}^{\dagger}. \tag*{(99)} ⟶ [ INNER † ⟶ reflect and square ] ⟶ OUTER † . (99) Equation (A4) maps the fixed root Majorana to γ u σ x \gamma_{u}\sigma_{x} γ u σ x O x † O x O_{x}^{\dagger}O_{x} O x † O x (72) .
For any chosen implementation convention, its Toffoli count can be recorded without exposing the two black boxes as
C B E = ν o u t C o u t + ν i n C i n + ∑ k = 0 K − 1 ν a n g , k C a n g , k ( J ) + C f i x e d ( T , N , J , b r o t ) , (100) C_{\mathrm{BE}} = \nu_{\mathrm{out}} C_{\mathrm{out}} + \nu_{\mathrm{in}} C_{\mathrm{in}} + \sum_{k=0}^{K-1} \nu_{\mathrm{ang},k} C_{\mathrm{ang},k}(J) + C_{\mathrm{fixed}}(\mathcal{T},N,J,b_{\mathrm{rot}}), \tag*{(100)} C BE = ν out C out + ν in C in + k = 0 ∑ K − 1 ν ang , k C ang , k ( J ) + C fixed ( T , N , J , b rot ) , (100) where the fixed integers ν o u t \nu_{\mathrm{out}} ν out ν i n \nu_{\mathrm{in}} ν in ν a n g , k \nu_{\mathrm{ang},k} ν ang , k A N G L E k \mathrm{ANGLE}_{k} ANGLE k 2 C a n g p a s s ( J ) 2C_{\mathrm{ang}}^{\mathrm{pass}}(J) 2 C ang pass ( J ) C a n g , k ( J ) C_{\mathrm{ang},k}(J) C ang , k ( J ) C o u t C_{\mathrm{out}} C out C i n C_{\mathrm{in}} C in
1/x \sqrt{x} x
Polynomial rescaling to cover x n o r m ∈ [ 2 , 4 ) x_{\mathrm{norm}} \in[2,4) x norm ∈ [ 2 , 4 )
The range of our interpolation is in [ 1 , 2 ) [1,2) [ 1 , 2 ) x n o r m x_{\mathrm{norm}} x norm [ 1 , 4 ) [1,4) [ 1 , 4 ) x n o r m ∈ [ 2 , 4 ) x_{\mathrm{norm}} \in[2,4) x norm ∈ [ 2 , 4 ) [ 1 , 2 ) [1,2) [ 1 , 2 ) u = x n o r m / 2 u = x_{\mathrm{norm}}/2 u = x norm /2
1 x n o r m = 1 2 u = 1 2 1 u (101) \frac{1}{\sqrt{x_{\mathrm{norm}}}} = \frac{1}{\sqrt{2u}} = \frac{1}{\sqrt{2}}\frac{1}{\sqrt{u}} \tag*{(101)} x norm 1 = 2 u 1 = 2 1 u 1 (101) and thus for x n o r m ∈ [ 2 , 4 ) x_{\mathrm{norm}} \in[2,4) x norm ∈ [ 2 , 4 ) (15) and Eq. (16) . To be more precise about the coefficients we need to use, note that for the upper ranges and by Eq. (101) we need to implement
1 x n o r m ≈ 1 2 P ( u ) (102) \frac{1}{\sqrt{x_{\mathrm{norm}}}} \approx\frac{1}{\sqrt{2}}P(u) \tag*{(102)} x norm 1 ≈ 2 1 P ( u ) (102) where P ( u ) P(u) P ( u ) (15) or Eq. (16) . To see how we need to modify the coefficients we can plug in our expression for u = x n o r m / 2 u = x_{\mathrm{norm}}/2 u = x norm /2
1 x n o r m = 1 2 [ a 0 − a 1 ( x n o r m 2 − 1 ) + a 2 ( x n o r m 2 − 1 ) 2 − a 3 ( x n o r m 2 − 1 ) 3 ] (103) \frac{1}{\sqrt{x_{\mathrm{norm}}}} = \frac{1}{\sqrt{2}}\left[a_{0} - a_{1}\left(\frac{x_{\mathrm{norm}}}{2}-1\right) + a_{2}\left(\frac{x_{\mathrm{norm}}}{2}-1\right)^{2} - a_{3}\left(\frac{x_{\mathrm{norm}}}{2}-1\right)^{3}\right] \tag*{(103)} x norm 1 = 2 1 [ a 0 − a 1 ( 2 x norm − 1 ) + a 2 ( 2 x norm − 1 ) 2 − a 3 ( 2 x norm − 1 ) 3 ] (103) = a 0 2 − a 1 2 2 ( x n o r m − 2 ) + a 2 4 2 ( x n o r m − 2 ) 2 − a 3 8 2 ( x n o r m − 2 ) 3 (104) = \frac{a_{0}}{\sqrt{2}} - \frac{a_{1}}{2\sqrt{2}}(x_{\mathrm{norm}}-2) + \frac{a_{2}}{4\sqrt{2}}(x_{\mathrm{norm}}-2)^{2} - \frac{a_{3}}{8\sqrt{2}}(x_{\mathrm{norm}}-2)^{3} \tag*{(104)} = 2 a 0 − 2 2 a 1 ( x norm − 2 ) + 4 2 a 2 ( x norm − 2 ) 2 − 8 2 a 3 ( x norm − 2 ) 3 (104) it is noted in Ref. [9 ] that when we scale the square root for some power of two the coefficients need to be scaled as
a 0 → a 0 / 2 m / 2 a_{0} \to a_{0}/2^{m/2} a 0 → a 0 / 2 m /2
a 1 → a 1 / 2 3 m / 2 a_{1} \to a_{1}/2^{3m/2} a 1 → a 1 / 2 3 m /2
a 2 → a 2 / 2 5 m / 2 a_{2} \to a_{2}/2^{5m/2} a 2 → a 2 / 2 5 m /2
a 3 → a 3 / 2 7 m / 2 a_{3} \to a_{3}/2^{7m/2} a 3 → a 3 / 2 7 m /2
with ( x − 1 ) (x-1) ( x − 1 ) ( x − 2 m ) (x-2^{m}) ( x − 2 m ) m = 1 m=1 m = 1
a 0 / 2 1 / 2 ⟹ a 0 / 2 a_{0}/2^{1/2} \Longrightarrow a_{0}/\sqrt{2} a 0 / 2 1/2 ⟹ a 0 / 2
a 1 / 2 3 / 2 ⟹ a 1 / 2 2 a_{1}/2^{3/2} \Longrightarrow a_{1}/2\sqrt{2} a 1 / 2 3/2 ⟹ a 1 /2 2
a 2 / 2 5 / 2 ⟹ a 2 / 4 2 a_{2}/2^{5/2} \Longrightarrow a_{2}/4\sqrt{2} a 2 / 2 5/2 ⟹ a 2 /4 2
a 3 / 2 7 / 2 ⟹ a 3 / 8 2 a_{3}/2^{7/2} \Longrightarrow a_{3}/8\sqrt{2} a 3 / 2 7/2 ⟹ a 3 /8 2
and ( x − 1 ) (x-1) ( x − 1 ) ( x n o r m − 2 ) (x_{\mathrm{norm}}-2) ( x norm − 2 ) b b b b b b 2 2 2 a a a x x x u ∈ [ 1.5 , 2 ] u \in[1.5,2] u ∈ [ 1.5 , 2 ]
P b ( u ) = b 0 − b 1 ( u − 1.5 ) + b 2 ( u − 1.5 ) 2 − b 3 ( u − 1.5 ) 3 . P_b(u)=b_0-b_1(u-1.5)+b_2(u-1.5)^2-b_3(u-1.5)^3. P b ( u ) = b 0 − b 1 ( u − 1.5 ) + b 2 ( u − 1.5 ) 2 − b 3 ( u − 1.5 ) 3 . We want to evaluate 1 / x 1/\sqrt{x} 1/ x x = u ⋅ 2 m x=u\cdot2^m x = u ⋅ 2 m u = x 2 m u=\frac{x}{2^m} u = 2 m x 1 x = 1 2 m ⋅ 1 u \frac{1}{\sqrt{x}}=\frac{1}{\sqrt{2^m}}\cdot\frac{1}{\sqrt{u}} x 1 = 2 m 1 ⋅ u 1
1 x ≈ 1 2 m / 2 [ b 0 − b 1 ( x 2 m − 1.5 ) + b 2 ( x 2 m − 1.5 ) 2 − b 3 ( x 2 m − 1.5 ) 3 ] . \frac{1}{\sqrt{x}}\approx\frac{1}{2^{m/2}}\left[b_0-b_1\left(\frac{x}{2^m}-1.5\right)+b_2\left(\frac{x}{2^m}-1.5\right)^2-b_3\left(\frac{x}{2^m}-1.5\right)^3\right]. x 1 ≈ 2 m /2 1 [ b 0 − b 1 ( 2 m x − 1.5 ) + b 2 ( 2 m x − 1.5 ) 2 − b 3 ( 2 m x − 1.5 ) 3 ] . Now, find a common denominator inside the parentheses. Notice that 1.5 1.5 1.5 1.5 ⋅ 2 m 2 m \frac{1.5\cdot2^m}{2^m} 2 m 1.5 ⋅ 2 m
( x 2 m − 1.5 ) = x − 1.5 ⋅ 2 m 2 m . \left(\frac{x}{2^m}-1.5\right)=\frac{x-1.5\cdot2^m}{2^m}. ( 2 m x − 1.5 ) = 2 m x − 1.5 ⋅ 2 m . Pulling those 2 m 2^m 2 m
1 x ≈ 1 2 m / 2 [ b 0 − b 1 2 m ( x − 1.5 ⋅ 2 m ) + b 2 2 2 m ( x − 1.5 ⋅ 2 m ) 2 − b 3 2 3 m ( x − 1.5 ⋅ 2 m ) 3 ] . \frac{1}{\sqrt{x}}\approx\frac{1}{2^{m/2}}\left[b_0-\frac{b_1}{2^m}(x-1.5\cdot2^m)+\frac{b_2}{2^{2m}}(x-1.5\cdot2^m)^2-\frac{b_3}{2^{3m}}(x-1.5\cdot2^m)^3\right]. x 1 ≈ 2 m /2 1 [ b 0 − 2 m b 1 ( x − 1.5 ⋅ 2 m ) + 2 2 m b 2 ( x − 1.5 ⋅ 2 m ) 2 − 2 3 m b 3 ( x − 1.5 ⋅ 2 m ) 3 ] . Finally, we distribute the 1 2 m / 2 \frac{1}{2^{m/2}} 2 m /2 1
1 x ≈ b 0 2 m / 2 − b 1 2 3 m / 2 ( x − 1.5 ⋅ 2 m ) + b 2 2 5 m / 2 ( x − 1.5 ⋅ 2 m ) 2 − b 3 2 7 m / 2 ( x − 1.5 ⋅ 2 m ) 3 . \frac{1}{\sqrt{x}}\approx\frac{b_0}{2^{m/2}}-\frac{b_1}{2^{3m/2}}(x-1.5\cdot2^m)+\frac{b_2}{2^{5m/2}}(x-1.5\cdot2^m)^2-\frac{b_3}{2^{7m/2}}(x-1.5\cdot2^m)^3. x 1 ≈ 2 m /2 b 0 − 2 3 m /2 b 1 ( x − 1.5 ⋅ 2 m ) + 2 5 m /2 b 2 ( x − 1.5 ⋅ 2 m ) 2 − 2 7 m /2 b 3 ( x − 1.5 ⋅ 2 m ) 3 . The scaling factors applied to b 0 , b 1 , b 2 , b 3 b_0,b_1,b_2,b_3 b 0 , b 1 , b 2 , b 3 a a a
b 0 → b 0 / 2 m / 2 b_0\rightarrow b_0/2^{m/2} b 0 → b 0 / 2 m /2
b 1 → b 1 / 2 3 m / 2 b_1\rightarrow b_1/2^{3m/2} b 1 → b 1 / 2 3 m /2
b 2 → b 2 / 2 5 m / 2 b_2\rightarrow b_2/2^{5m/2} b 2 → b 2 / 2 5 m /2
b 3 → b 3 / 2 7 m / 2 b_3\rightarrow b_3/2^{7m/2} b 3 → b 3 / 2 7 m /2
And ( x − 1.5 ) (x-1.5) ( x − 1.5 ) ( x − 1.5 ⋅ 2 m ) (x-1.5\cdot2^m) ( x − 1.5 ⋅ 2 m )
Example calculation of ∣ s ⟩ \lvert s\rangle ∣ s ⟩
We illustrate the normalization stage on the n p = 6 n_p=6 n p = 6 x 13 … x 0 x_{13}\ldots x_0 x 13 … x 0 n p n_p n p
Starting with the highest two bits ( x 13 , x 12 ) (x_{13},x_{12}) ( x 13 , x 12 ) f 0 f_0 f 0 n p n_p n p
Controlled on f 0 f_0 f 0 ( x 0 , x 2 , … , x 12 ) (x_0,x_2,\ldots,x_{12}) ( x 0 , x 2 , … , x 12 ) ( x 1 , x 3 , … , x 13 ) (x_1,x_3,\ldots,x_{13}) ( x 1 , x 3 , … , x 13 ) L L L L − 1 L-1 L − 1
Why cyclic rotation leaves clean zeros: Whenever f 0 = 1 f_0=1 f 0 = 1 ( x 13 , x 12 ) (x_{13},x_{12}) ( x 13 , x 12 ) x 11 → x 13 , … , x 1 → x 3 x_{11}\to x_{13},\ldots,x_1\to x_3 x 11 → x 13 , … , x 1 → x 3 x 13 x_{13} x 13 x 1 x_1 x 1 x 12 x_{12} x 12 x 0 x_0 x 0
Cascading the shifts (k > 0 k>0 k > 0 We repeat this process for k ∈ { 1 , … , n p − 1 } k\in\{1,\ldots,n_p-1\} k ∈ { 1 , … , n p − 1 } k k k ( x 2 n p + 1 , x 2 n p ) (x_{2n_p+1},x_{2n_p}) ( x 2 n p + 1 , x 2 n p ) f k f_k f k 2 k 2k 2 k
5. Truncated Loop Optimization: Because a shift at step k k k < k < k < k 2 k 2k 2 k x x x k k k n p + 1 − k n_p + 1 - k n p + 1 − k n p − k n_p - k n p − k n p = 6 n_p = 6 n p = 6
Step k = 0 k = 0 k = 0 f 0 f_0 f 0 x 13 … x 0 x_{13}\ldots x_0 x 13 … x 0 → \to →
Step k = 1 k = 1 k = 1 f 1 f_1 f 1 x 13 … x 2 x_{13}\ldots x_2 x 13 … x 2 → \to →
Step k = 2 k = 2 k = 2 f 2 f_2 f 2 x 13 … x 4 x_{13}\ldots x_4 x 13 … x 4 → \to →
Step k = 3 k = 3 k = 3 f 3 f_3 f 3 x 13 … x 6 x_{13}\ldots x_6 x 13 … x 6 → \to →
Step k = 4 k = 4 k = 4 f 4 f_4 f 4 x 13 … x 8 x_{13}\ldots x_8 x 13 … x 8 → \to →
Step k = 5 k = 5 k = 5 f 5 f_5 f 5 x 13 … x 10 x_{13}\ldots x_{10} x 13 … x 10 → \to →
Summing over k = 0 , … , n p − 1 k = 0,\ldots,n_p - 1 k = 0 , … , n p − 1 n p ( n p + 1 ) n_p(n_p + 1) n p ( n p + 1 ) n p = 6 n_p = 6 n p = 6 n p = 10 n_p = 10 n p = 10
6. Total Complexity: Adding the n p n_p n p f k f_k f k n p ( n p + 2 ) n_p(n_p + 2) n p ( n p + 2 ) n p = 6 n_p = 6 n p = 6 n p ( n p + 1 ) n_p(n_p + 1) n p ( n p + 1 ) n p n_p n p x n o r m x_{\mathrm{norm}} x norm n p n_p n p ∣ f ⟩ \lvert f\rangle ∣ f ⟩ s s s
Details on cubic polynomial approximation for 1 / w 1/\sqrt{w} 1/ w
The different polynomial coefficients are input data that must be output coherently via QROM. The selection state for QROM will be the three most significant bits of x n o r m x_{\mathrm{norm}} x norm x n o r m = b 13 b 12 . b 11 b 10 … x_{\mathrm{norm}} = b_{13}b_{12}.b_{11}b_{10}\ldots x norm = b 13 b 12 . b 11 b 10 … b 13 = 1 b_{13} = 1 b 13 = 1 x n o r m ∈ [ 2 , 4 ) x_{\mathrm{norm}} \in[2,4) x norm ∈ [ 2 , 4 )
b 13 = 1 b_{13} = 1 b 13 = 1 b 12 = 0 b_{12} = 0 b 12 = 0 x n o r m ∈ [ 2 , 3 ) x_{\mathrm{norm}} \in[2,3) x norm ∈ [ 2 , 3 ) (15) scaled by something proportional to 1 / 2 1/\sqrt{2} 1/ 2
b 13 = 1 b_{13} = 1 b 13 = 1 b 12 = 1 b_{12} = 1 b 12 = 1 x n o r m ∈ [ 3 , 4 ) x_{\mathrm{norm}} \in[3,4) x norm ∈ [ 3 , 4 ) (16) scaled by something proportional to 1 / 2 1/\sqrt{2} 1/ 2
b 13 = 0 b_{13} = 0 b 13 = 0 b 12 = 1 b_{12} = 1 b 12 = 1 b 11 = 0 b_{11} = 0 b 11 = 0 x n o r m ∈ [ 1 , 1.5 ) x_{\mathrm{norm}} \in[1,1.5) x norm ∈ [ 1 , 1.5 ) (15)
b 13 = 0 b_{13} = 0 b 13 = 0 b 11 = 1 b_{11} = 1 b 11 = 1 x n o r m ∈ [ 1.5 , 2 ) x_{\mathrm{norm}} \in[1.5,2) x norm ∈ [ 1.5 , 2 ) (16)
The QROM cost for these four entries is 2 Toffolis because it is not controlled (two AND gates decode t a = ¬ b 13 ∧ b 11 t_a = \neg b_{13} \land b_{11} t a = ¬ b 13 ∧ b 11 t b = b 13 ∧ b 12 t_b = b_{13} \land b_{12} t b = b 13 ∧ b 12 { 1 , b 13 , t a , t b } \{1,b_{13},t_a,t_b\} { 1 , b 13 , t a , t b } ∣ c 0 ⟩ , ∣ c 1 ⟩ , ∣ c 2 ⟩ , ∣ c 3 ⟩ \lvert c_0\rangle,\lvert c_1\rangle,\lvert c_2\rangle,\lvert c_3\rangle ∣ c 0 ⟩ , ∣ c 1 ⟩ , ∣ c 2 ⟩ , ∣ c 3 ⟩
Next, we describe a series of optimizations for evaluating the polynomial once the coefficients are obtained by the QROM read. We need to evaluate these polynomials on x n o r m x_{\mathrm{norm}} x norm (15) and Eq. (16) it would seem we need many multiplications and subtractions. The computation can be significantly simplified by rewriting it with Horner’s rule as (for the case of the range [ 1 , 3 / 2 ) [1,3/2) [ 1 , 3/2 )
a 0 − ( x − 1 ) { a 1 − ( x − 1 ) [ a 2 − a 3 ( x − 1 ) ] } = a 0 − a 1 ( x − 1 ) + a 2 ( x − 1 ) 2 − a 3 ( x − 1 ) 3 . (105) a_0-(x-1)\{a_1-(x-1)[a_2-a_3(x-1)]\}=a_0-a_1(x-1)+a_2(x-1)^2-a_3(x-1)^3. \tag*{(105)} a 0 − ( x − 1 ) { a 1 − ( x − 1 ) [ a 2 − a 3 ( x − 1 )]} = a 0 − a 1 ( x − 1 ) + a 2 ( x − 1 ) 2 − a 3 ( x − 1 ) 3 . (105) and for the range [ 3 / 2 , 2 ) [3/2,2) [ 3/2 , 2 )
b 0 − ( x − 1.5 ) { b 1 − ( x − 1.5 ) [ b 2 − b 3 ( x − 1.5 ) ] } = b 0 − b 1 ( x − 1.5 ) + b 2 ( x − 1.5 ) 2 − b 3 ( x − 1.5 ) 3 . (106) b_0-(x-1.5)\{b_1-(x-1.5)[b_2-b_3(x-1.5)]\}=b_0-b_1(x-1.5)+b_2(x-1.5)^2-b_3(x-1.5)^3. \tag*{(106)} b 0 − ( x − 1.5 ) { b 1 − ( x − 1.5 ) [ b 2 − b 3 ( x − 1.5 )]} = b 0 − b 1 ( x − 1.5 ) + b 2 ( x − 1.5 ) 2 − b 3 ( x − 1.5 ) 3 . (106) The above rewrites of the equation need Δ x = x − 1 \Delta x=x-1 Δ x = x − 1 Δ x = x − 1.5 \Delta x=x-1.5 Δ x = x − 1.5 Δ x = x − 2 \Delta x=x-2 Δ x = x − 2 Δ x = x − 3 \Delta x=x-3 Δ x = x − 3 P a ( u ) P_a(u) P a ( u ) P b ( u ) P_b(u) P b ( u )
case 3: x n o r m ∈ [ 1 , 1.5 ) x_{\mathrm{norm}}\in[1,1.5) x norm ∈ [ 1 , 1.5 ) x n o r m = u x_{\mathrm{norm}}=u x norm = u Δ x = ( x − 1 ) \Delta x=(x-1) Δ x = ( x − 1 ) x n o r m x_{\mathrm{norm}} x norm 01.0 b 10 b 9 … 01.0b_{10}b_9\ldots 01.0 b 10 b 9 … b 13 = 0 , b 12 = 1 , b 11 = 0 b_{13}=0,b_{12}=1,b_{11}=0 b 13 = 0 , b 12 = 1 , b 11 = 0 b 12 b_{12} b 12 Δ x \Delta x Δ x 00.0 b 10 b 9 … 00.0b_{10}b_9\ldots 00.0 b 10 b 9 …
case 4: x n o r m ∈ [ 1.5 , 2 ) x_{\mathrm{norm}}\in[1.5,2) x norm ∈ [ 1.5 , 2 ) x n o r m = u x_{\mathrm{norm}}=u x norm = u Δ x = ( x − 1.5 ) \Delta x=(x-1.5) Δ x = ( x − 1.5 ) x n o r m x_{\mathrm{norm}} x norm 01.1 b 10 b 9 … 01.1b_{10}b_9\ldots 01.1 b 10 b 9 … b 13 = 0 , b 12 = 1 , b 11 = 1 b_{13}=0,b_{12}=1,b_{11}=1 b 13 = 0 , b 12 = 1 , b 11 = 1 b 12 b_{12} b 12 b 11 b_{11} b 11 Δ x \Delta x Δ x 00.0 b 10 b 9 … 00.0b_{10}b_9\ldots 00.0 b 10 b 9 …
case 1: x n o r m ∈ [ 2 , 3 ) x_{\mathrm{norm}} \in[2,3) x norm ∈ [ 2 , 3 ) x n o r m x_{\mathrm{norm}} x norm u = x n o r m / 2 u=x_{\mathrm{norm}}/2 u = x norm /2 u ∈ [ 1 , 1.5 ) u \in[1,1.5) u ∈ [ 1 , 1.5 ) P a ( u ) P_a(u) P a ( u )
( u − 1 ) = ( x n o r m 2 − 1 ) = 1 2 ( x n o r m − 2 ) (107) (u-1)=\left(\frac{x_{\mathrm{norm}}}{2}-1\right)=\frac{1}{2}(x_{\mathrm{norm}}-2) \tag*{(107)} ( u − 1 ) = ( 2 x norm − 1 ) = 2 1 ( x norm − 2 ) (107) which is the same calculation as Eq. (103) . As you can see each of the powers of 2 gets absorbed into a redefinition of the coefficient a i a_i a i x n o r m x_{\mathrm{norm}} x norm 10. b 11 b 10 … 10.b_{11}b_{10}\ldots 10. b 11 b 10 … b 13 = 1 b_{13}=1 b 13 = 1 b 12 = 0 b_{12}=0 b 12 = 0 b 13 b_{13} b 13 Δ x \Delta x Δ x 00. b 11 b 10 … 00.b_{11}b_{10}\ldots 00. b 11 b 10 …
case 2: x n o r m ∈ [ 3 , 4 ) x_{\mathrm{norm}} \in[3,4) x norm ∈ [ 3 , 4 ) x n o r m x_{\mathrm{norm}} x norm u = x n o r m / 2 u=x_{\mathrm{norm}}/2 u = x norm /2 u ∈ [ 1.5 , 2 ) u \in[1.5,2) u ∈ [ 1.5 , 2 )
u − 1.5 = ( x n o r m 2 − 1.5 ) = 1 2 ( x n o r m − 3 ) (108) u-1.5=\left(\frac{x_{\mathrm{norm}}}{2}-1.5\right)=\frac{1}{2}(x_{\mathrm{norm}}-3) \tag*{(108)} u − 1.5 = ( 2 x norm − 1.5 ) = 2 1 ( x norm − 3 ) (108) which matches the scaling coefficients and Δ x \Delta x Δ x P b ( u ) P_b(u) P b ( u ) x n o r m x_{\mathrm{norm}} x norm 11. b 11 b 10 … 11.b_{11}b_{10}\ldots 11. b 11 b 10 … b 13 = 1 b_{13}=1 b 13 = 1 b 12 = 1 b_{12}=1 b 12 = 1 b 13 b_{13} b 13 b 12 b_{12} b 12 Δ x \Delta x Δ x 00. b 11 b 10 … 00.b_{11}b_{10}\ldots 00. b 11 b 10 …
We need the bits of Δ x \Delta x Δ x x n o r m x_{\mathrm{norm}} x norm b 13 b_{13} b 13 b 12 b_{12} b 12 Δ x \Delta x Δ x b 10 b_{10} b 10 b 0 b_0 b 0 x n o r m x_{\mathrm{norm}} x norm b 11 b_{11} b 11 b 11 b_{11} b 11 Δ x \Delta x Δ x
Case 1 ([2, 3)): We subtract 2. We keep the 0.5s place. ⟹ \Longrightarrow ⟹ b 11 b_{11} b 11
Case 2 ([3, 4)): We subtract 3. We keep the 0.5s place. ⟹ \Longrightarrow ⟹ b 11 b_{11} b 11
Case 3 ([1, 1.5)): We subtract 1. b 11 b_{11} b 11 Δ x 11 = 0 \Delta x_{11}=0 Δ x 11 = 0
Case 4 ([1.5, 2)): We subtract 1.5, clearing the 0.5s place. Even though b 11 = 1 b_{11}=1 b 11 = 1 x n o r m x_{\mathrm{norm}} x norm Δ x 11 \Delta x_{11} Δ x 11
Look at the logic for b 11 b_{11} b 11 b 11 b_{11} b 11 b 13 = 1 b_{13}=1 b 13 = 1 b 13 = 0 b_{13}=0 b 13 = 0 Δ x 11 = b 13 ∧ b 11 \Delta x_{11}=b_{13}\land b_{11} Δ x 11 = b 13 ∧ b 11 ∣ Δ x ⟩ \lvert\Delta x\rangle ∣ Δ x ⟩ b 13 ∧ b 11 b_{13}\land b_{11} b 13 ∧ b 11 Δ x \Delta x Δ x x n o r m x_{\mathrm{norm}} x norm
Finally, we perform the Horner multiplications and the Newton-Raphson step. The polynomial to evaluate is
y 0 = a 0 − a 1 Δ x + a 2 Δ x 2 − a 3 Δ x 3 = a 0 − Δ x ⋅ ( a 1 − Δ x ⋅ ( a 2 − a 3 ⋅ Δ x ) ) , (109) y_0=a_0-a_1\Delta x+a_2\Delta x^2-a_3\Delta x^3=a_0-\Delta x\cdot\left(a_1-\Delta x\cdot(a_2-a_3\cdot\Delta x)\right), \tag*{(109)} y 0 = a 0 − a 1 Δ x + a 2 Δ x 2 − a 3 Δ x 3 = a 0 − Δ x ⋅ ( a 1 − Δ x ⋅ ( a 2 − a 3 ⋅ Δ x ) ) , (109) which requires three multiplications and three subtractions. In our reference circuit, we use 14 bits of Δ x \Delta x Δ x n p = 6 n_p=6 n p = 6 n p = 10 n_p=10 n p = 10 ∣ c 0 ⟩ , ∣ c 1 ⟩ , ∣ c 2 ⟩ , ∣ c 3 ⟩ \lvert c_0\rangle,\lvert c_1\rangle,\lvert c_2\rangle,\lvert c_3\rangle ∣ c 0 ⟩ , ∣ c 1 ⟩ , ∣ c 2 ⟩ , ∣ c 3 ⟩ y 0 y_0 y 0 z = x n o r m y 0 z=x_{\mathrm{norm}}y_0 z = x norm y 0 3 − z 3-z 3 − z p = y 0 ( 3 − z ) p=y_0(3-z) p = y 0 ( 3 − z ) 30 + n p 30+n_p 30 + n p n p = 10 n_p=10 n p = 10 n p = 6 n_p=6 n p = 6 y 1 = p ≪ s y_1=p\ll s y 1 = p ≪ s
In-place Horner evaluation: Rather than allocating six fresh output registers, the circuit updates the coefficient registers in place over three steps: c 2 = T ( c 3 Δ x ) c_2=T(c_3\Delta x) c 2 = T ( c 3 Δ x ) c 1 = T ( c 2 Δ x ) c_1=T(c_2\Delta x) c 1 = T ( c 2 Δ x ) c 0 = T ( c 1 Δ x ) c_0=T(c_1\Delta x) c 0 = T ( c 1 Δ x ) y 0 y_0 y 0 ∣ c 0 ⟩ \lvert c_0\rangle ∣ c 0 ⟩ n p = 10 n_p=10 n p = 10 n p = 6 n_p=6 n p = 6 T ( ⋅ ) T(\cdot) T ( ⋅ ) Δ x \Delta x Δ x
Mantissa multiply and correction: Multiplies T ( y 0 2 ) T(y_0^2) T ( y 0 2 ) 2 n p + 2 2n_p+2 2 n p + 2 x n o r m x_{\mathrm{norm}} x norm z = T ( x n o r m ⋅ T ( y 0 2 ) ) z=T(x_{\mathrm{norm}}\cdot T(y_0^2)) z = T ( x norm ⋅ T ( y 0 2 )) n p = 10 n_p=10 n p = 10 n p = 6 n_p=6 n p = 6 z z z 3 − z 3-z 3 − z X X X ( 3 ≪ 28 ) + 1 (3\ll28)+1 ( 3 ≪ 28 ) + 1
Final NR multiply and denormalization: Computes p = T ( y 0 ⋅ z ) ≈ 2 / x n o r m p=T(y_0\cdot z)\approx2/\sqrt{x_{\mathrm{norm}}} p = T ( y 0 ⋅ z ) ≈ 2/ x norm y 1 y_1 y 1 y 1 = p ≪ s y_1=p\ll s y 1 = p ≪ s y 1 [ k : k + 31 ] y_1[k:k+31] y 1 [ k : k + 31 ] f k f_k f k k = 0 , … , n p − 1 k=0,\ldots,n_p-1 k = 0 , … , n p − 1 30 n p 30n_p 30 n p n p = 10 n_p=10 n p = 10 n p = 6 n_p=6 n p = 6
Phasing and uncomputation: Applies one Z_POW gate of 2 j − F 2^{j-F} 2 j − F j ∈ { 0 , … , F − 1 } j \in\{0, \ldots, F - 1\} j ∈ { 0 , … , F − 1 } y 1 y_1 y 1 F = 30 + n p F = 30 + n_p F = 30 + n p ∣ 0 ⟩ |0\rangle ∣0 ⟩ ∣ x ⟩ |x\rangle ∣ x ⟩
