Simulation

qirc program.ll compiles the program and runs it. Which simulator runs it depends on the program: small programs use a state vector, and wider ones use whichever structured kernel fits the gates they contain.

The first row that fits runs the program, and the run prints which kernel it used.

State vector

The state vector is stored as separate real and imaginary arrays, so a single qubit gate is a 2x2 complex matrix applied to every pair of amplitudes at once. On x86_64 with AVX2 and FMA, four amplitudes are processed per instruction.

When the target qubit is 2 or higher, each pair’s two halves are contiguous and load directly. When the target is qubit 0 or 1, both halves sit in the same register and are paired with a permute. Controls are a bit mask: a control above bit 1 is constant across a register, so whole registers are skipped, and a control on bit 0 or 1 is blended per lane. This covers Toffoli and controlled swap without a separate code path. CPUs without AVX2 use a scalar fallback. From 21 qubits each gate is split across the CPU’s threads, a chunk of the amplitudes each, or a pair of chunks when the target qubit is above the chunk size. At 26 qubits, where the state is a gigabyte and every gate reads and writes all of it, 16 threads make a random circuit 1.5 times faster, which is as far as memory bandwidth allows without applying several gates in one pass.

Straight line programs are evolved once and sampled. Programs that branch on a measurement, reset a qubit, or act on a qubit after measuring it are simulated shot by shot with real collapse.

Wide programs

A program with more than 20 qubits whose gates are all Clifford (H, S, the Paulis, SX, CNOT, CZ, CY, swap, and rotations by multiples of a quarter turn) runs on a stabilizer tableau instead, up to 5,000 qubits, with branches, resets and measurements. The tableau is stored row by row as bit words, so combining two rows costs a few popcounts, and a straight line program is prepared once and each shot measures a copy. A 1,000 qubit GHZ state takes under a second for 1,000 shots. qirc diff uses the same tableau for such programs, following each measurement that is not already determined and comparing final states through the reduced row echelon form of their stabilizer groups, so it can check a compile of a large Clifford circuit exactly.

A program with more than 20 qubits whose gates only permute basis states (X, Y, Toffoli and other multi controlled X gates, swaps) or only add phases (Z, S, T, Rz, CZ and other controlled phases) runs on a plain bit vector instead, with no limit on its width, since such gates keep a basis state a basis state. Arithmetic like a ripple carry adder falls in this class even though its Toffolis are not Clifford: two 5000 bit numbers add on 10,002 qubits in 64 ms, and noise still applies, as a bit flip for X and Y errors and nothing for Z errors, which only change the phase. The run prints the circuit only when it has at most 200,000 cells of qubits by depth, and points to --emit circuit otherwise.

Cutting also samples. A straight line program above 30 qubits that measures at most 20 of them at the end, and cuts within the same limits, runs as circuit cutting before falling back to the matrix product state: each piece’s states give a table of overlaps for every outcome of its measured qubits, and the tables recombine into the exact distribution over the measured bits, which the shots are drawn from. The 40 qubit circuit above, measuring 6 qubits, samples 4,000 shots in about 2 seconds.

A straight line program with more than 30 qubits whose gates are Clifford apart from up to 10 rotations by other angles, such as T gates, runs on a sum of stabilizer states. Each state is kept in CH form, which unlike a tableau tracks its global phase, so the terms of the sum can interfere. Lowered to rz, sx and cx, a CNOT or a quarter turn is a row operation on each term, a Hadamard folds the two basis strings it produces back into one state with a few CNOTs, CZs and an S, and a rotation by any other angle splits every term into two with cos(θ/2) I - i sin(θ/2) Z. Shots are drawn exactly by gate by gate sampling, which needs only amplitudes: every shot moves through the circuit with the state, a CNOT or X flips its bits, a phase changes nothing, and only an sx redraws one bit from the amplitudes of its two neighbours. Shots that agree so far share those amplitudes. Measurements and resets in the middle of a straight line program are deferred: a measurement becomes a CNOT onto a fresh qubit read at the end, and a reset swaps the qubit with a fresh one, which gives the same statistics exactly. A program that branches on a result goes to the matrix product state instead. A 40 qubit circuit with three T gates spread along a CNOT chain takes 3 seconds for 40,000 shots, and every qubit’s frequency of 1 matches the exact value from the light cone expectation to within 0.0004.

Any other program with more than 30 qubits runs as a matrix product state: one small tensor per qubit, joined by bonds whose size grows with entanglement. A two qubit gate on neighbouring tensors is applied to their product and split back with a singular value decomposition, gates between distant qubits move one tensor next to the other with swaps, and the bonds are capped at 32 by default or at --bond n, dropping the smallest singular values. The run reports the weight it dropped, so an approximate answer says so. Without the cap the result is exact, which the tests check against the state vector amplitude by amplitude. A 60 qubit circuit of four layers of arbitrary rotations and CNOTs runs in 1.3 seconds and its marginals match exact light cone values to within sampling noise, and a 40 qubit circuit of twelve layers drops 0.09% of the weight at bond 32. --bond also forces this kernel on a small program.

Expectation values

--observable prints the exact expectation value of a sum of Pauli strings, written like Z0 Z1 + 0.5 X2 - 2*Y3. Measurements at the end of the program are left out, so the value is taken on the state just before them, and a program that branches on mid circuit measurements is averaged exactly over every branch, weighted by its probability. On the teleportation example the teleported qubit gives 0.707107 for both Z2 and X2, as ry(pi/4) should. A program without mid circuit measurements is not limited to 30 qubits: each Pauli string only depends on the gates in its backward light cone, so each is simulated on a small circuit over just the qubits of its cone. A 200 qubit layer of ry rotations and CNOT pairs gives Z199 + X0 X1 exactly in 50 ms. A cone of permutation and phase gates is a basis state, where a string with an X or Y gives 0 and a string of Z gives the sign of the parity of its bits. A cone of Clifford gates is simulated on the stabilizer tableau instead, however wide it is: after basis changes and CNOTs fold the Pauli string onto one qubit, its value is 0 if measuring that qubit would be random and the sign of the deterministic outcome otherwise, so Y0 Y1 X2 ... X999 on a 1000 qubit GHZ state gives -1. A cone that is neither Clifford nor small enough for a state vector is lowered to rz, sx and cx and the Pauli string is pushed backwards through it in the Heisenberg picture: a CNOT maps each string to one string, a rotation splits a string that anticommutes with it into a cosine and a sine part, and at the start only strings of I and Z count. This stays exact while few rotations touch the cone, so rz(0.3), t and rx(0.2) on a 100 qubit GHZ state give X0 ... X99 + Z0 Z99 = 1.446627 in 80 ms, and it gives up once a term expands past 65536 strings.

When the expansion gives up, the cone is cut instead. qirc splits its qubits into the fewest groups of at most 30 that as few two qubit gates as possible connect, peeling off one group at a time and moving qubits between groups while that removes connecting gates, and tries smaller groups when larger ones would need more memory than it allows, and writes each connecting gate as a sum of products of one qubit operators: a controlled gate as |0⟩⟨0| ⊗ I + |1⟩⟨1| ⊗ U, and a swap as half the sum of P ⊗ P over the four Paulis. Each half is simulated as a state vector once for every choice of those terms, branching at each cut so the gates before it run only once, and the expectation is recombined exactly from the overlaps between those states, one for each pair of choices, so the work grows with the square of the number of choices and a few cuts are practical. A 40 qubit circuit made of two 20 qubit halves of random rotations and CNOT chains, joined by two CNOTs, gives Z0 Z39 + X10 + Z19 Z20 exactly in about 5 seconds, and followed by its own inverse, four cuts, it gives exactly 2. Terms whose light cones keep the same gates share one set of simulations. Three 20 qubit blocks joined by two CNOTs give an exact expectation on 60 qubits in about a second.