BQN's Quantum Noise

Preamble

We will implement and test a Schrödinger-style¹ quantum simulator in the BQN² programming language. Initially, we import the necessary system functions and define a 1-modifier for handling complex-valued matrix products. Next, we define a namespace containing various quantum gates:

Sin‿Cos‿GCD ← •math
U ← •rand.Range
_cp ← {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽}

g ← {
  IM ← (⊢⋈≢⥊⟜0⊢)¨
  x‿id‿h ⇐ (⌽‿⊢{𝕎𝕩}¨<=⌜˜↕2) ∾○IM <-⌾(1‿1⊸⊑)2‿2⥊÷√2
  swap‿cnot ⇐ IM ⟨1‿2, 2‿3⟩ {⌽⌾(𝕨⊸⊏)𝕩}¨ <=⌜˜↕4
  P ⇐ {⟨=⌜˜↕4,  4‿4⥊0⟩ {𝕨⌾(3‿3⊸⊑) 𝕩}¨˜ Sin⊸⋈⟜Cos 𝕩}
}

Interpreter

The (400 chars³) quantum interpreter is based on references arXiv:1711.02086 and arXiv:1608.03355. For simplicity, we always measure at the end of the execution:

Q ← {𝕊qb‿sc‿r:
  wf ← (1⌾⊑⋈⊢)⥊⟜0 2⋆qb
  M‿K ← ⟨+˝∘×⎉1‿∞ _cp, {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}⟩
  E ← {0𝕊𝕩:1; K⍟(𝕨-1)˜𝕩}
  L ← {K´ ⟨(qb-𝕨+⌊2⋆⁼≠⊑𝕩) E g.id, 𝕩, 𝕨 E g.id⟩}
  T ← ∾↕∘≠{a←𝕩, i←𝕨{𝕩⊑a}•_while_{𝕩<𝕨}𝕨⊑a, 𝕨<◶⟨⟩‿{(⊢∾1⊸↓∘⌽)𝕨+↕𝕩-𝕨}i}¨<
  A ← {qs‿gs𝕊v:
    1⊸=◶{𝕊𝕩: ui ← 0 L gs
      v M˜ {𝕊⟨⟩:ui; (M˜´M⟜ui⟜M˜M´) L⟜g.swap¨ 𝕩} T qs (⌽∘⊢∾¬∘∊/⊣)˜ ↕qb
    }‿(v M˜ gs L˜ ⊑qs) ≠qs}
  »⊸<∨` 0>r-`>+○(×˜)˝ wf A´ ⌽sc
}

Physical context

The interpreter’s logic is straightforward: it evolves the wavefunction \(\Psi_0^n\) in the full Hilbert space, whose dimension is \(2^n\) for \(n\) qubits. Each quantum gate is extended to this space through Kronecker lifting (E), performed by the function L:

\begin{equation*} L_U = I \otimes \cdots U \cdots \otimes I \end{equation*}

Note that the function L performs more than a plain Kronecker embedding. Before constructing \(L_U\), the affected qubits are permuted to become index-adjacent in the tensor order. This permutation is decomposed into adjacent swaps, applied as unitary transpositions (T). The system then evolves under the sequence of lifted gates, defined by A:

\begin{equation*} \Psi_m^n = \Psi_0^n \prod_i^m L_{U_i} \end{equation*}

Measurement is implemented by sampling from the CDF of the squared amplitudes of the state vector. Statistically, this corresponds to drawing a discrete random variable \(q_j\) whose cumulative distribution is:

\begin{equation*} F(q_j) = \sum_{i \le j}p(q_i) \end{equation*}

Given a uniform random variable \(U \sim \text{Uniform}(0, 1)\), we sample by inversion:

\begin{equation*} \text{Pr}(U \in (F(q_{j-1}), F(q_{j})]) = \text{Pr}(F(q_{j-1}) < U \le F(q_{j})) = F(q_{j}) - F(q_{j-1}) = p(q_{j}) \end{equation*}

Shor's algorithm

As a test case, we employ the quantum circuit of Shor's algorithm for the number fifteen and base eleven, following references arXiv:1804.03719 and arXiv:2306.09122. The resulting compiled circuit uses five qubits, three of which serve as control. To enhance statistical accuracy, the experiment is repeated multiple times.

n‿a‿qb‿r ← ⟨15, 11, 5, 0 U˜ 2⋆3⟩

sc ← ⟨
  ⟨0⟩‿g.h ⋄ ⟨1⟩‿g.h ⋄ ⟨2⟩‿g.h
  ⟨2, 3⟩‿g.cnot ⋄ ⟨2, 4⟩‿g.cnot ⋄ ⟨1⟩‿g.h
  ⟨⟨1, 0⟩, g.P π÷2⟩ ⋄ ⟨0⟩‿g.h
  ⟨⟨1, 2⟩, g.P π÷4⟩ ⋄ ⟨⟨0, 2⟩, g.P π÷2⟩ ⋄ ⟨2⟩‿g.h
⟩

C ← {n (⊣≡×´∘GCD) +‿-{𝕩𝕎1}¨ <a⋆(≠÷2×⊑∘⍒) 0⌾⊑+˝∘‿(2⋆qb-2)⥊𝕩}

The function C performs the classical post-processing part of the algorithm. Given the measurement registers, it reconstructs the candidate period \(r\) and checks the Shor certificate: it is true iff both \(\gcd(a^{\,r/2}\pm 1,N)\) are non-trivial and their product equals \(N\).

Let’s run the quantum simulation and complete the reasoning in the classical domain; let us kill the ignorabimus. Wir müssen wissen, wir werden wissen!⁴

C >+˝{Q qb‿sc‿𝕩}¨ r

Compare the result with that obtained using an IBM's Eagle QPU.

Epilogue

This was, of course, a conceptually simple simulation, the minimal instance of Shor’s algorithm that has been verified experimentally. The interpreter itself remains in its basic, non-optimized form. Yet for me, the goal was different: to stress-test BQN on a genuinely technical task. I enjoyed distilling the essence of a long and intricate implementation into a few concise expressions, like trimming a bonsai.

In that process, the language revealed its character. Why the hieroglyphs, you may ask? The tacit and functional style woven from combinators, turns programming into a quiet algebraic game. The resulting program feels less written than discovered, its complexity compressed into a handful of primitives. BQN is the epitome of minimalism's power, and the primitive statistics below tell the same story in numbers:

Primitive's stats

The prog string contains the full source code. We used:

prog (+´⊸≍⟜≠∊)˜ ⊑¨•primitives

⟨ 44 64 ⟩

With this distribution:

⍉>(⍷∾≠)¨∘(⊐⊸⊔∊/⊣)⟜(⊑¨•primitives)˜ prog

┌─                                                                                                                                                                                 
╵ '-' '´' '¨' '⋈' '+' '⟜' '⌽' '⊢' '≢' '⥊' '<' '=' '⌜' '˜' '↕' '∾' '○' '⌾' '⊸' '⊑' '÷' '√' '⊏' '⋆' '˝' '∘' '×' '⎉' '≡' '⊣' '⌊' '⁼' '≠' '⍟' '◶' '↓' '¬' '∊' '/' '»' '∨' '`' '>' '⍒'  
  8   8   10  5   8   3   6   7   1   5   9   6   3   12  6   5   2   5   7   9   5   1   1   5   4   8   5   1   3   3   1   1   5   1   2   1   1   1   1   1   1   2   3   1    
                                                                                                                                                                                  ┘

BQN is also fast:

Benchmarks

While the interpreter's performance is not particularly optimized, here is a comparison with the equivalent Common Lisp code:

Benchmark 1: cbqn -f ./bqn/q.bqn
  Time (mean ± σ):      5.468 s ±  0.077 s    [User: 5.427 s, System: 0.005 s]
  Range (min … max):    5.358 s …  5.535 s    5 runs
 
Benchmark 2: sbcl --script ../../supp/perf_qi/q.lisp
  Time (mean ± σ):     37.114 s ±  0.893 s    [User: 37.544 s, System: 0.207 s]
  Range (min … max):   36.457 s … 38.634 s    5 runs
 
Summary
  cbqn -f ./bqn/q.bqn ran
    6.79 ± 0.19 times faster than sbcl --script ../../supp/perf_qi/q.lisp

And here is a full program's profile. All time is spent in the Kronecker and matrix products:

)profile C >+˝{Q qb‿sc‿𝕩}¨ r

Got 25361 samples
(REPL): 25361 samples:
     2│  Q ← {𝕊qb‿sc‿r:
     1│    wf ← (1⌾⊑⋈⊢)⥊⟜0 2⋆qb
  2471│    M‿K ← ⟨+˝∘×⎉1‿∞ _cp, {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}⟩
    26│    E ← {0𝕊𝕩:1; K⍟(𝕨-1)˜𝕩}
    39│    L ← {K´ ⟨(qb-𝕨+⌊2⋆⁼≠⊑𝕩) E g.id, 𝕩, 𝕨 E g.id⟩}
    16│    T ← ∾↕∘≠{a←𝕩, i←𝕨{𝕩⊑a}•_while_{𝕩<𝕨}𝕨⊑a, 𝕨<◶⟨⟩‿{(⊢∾1⊸↓∘⌽)𝕨+↕𝕩-𝕨}i}¨<
     1│    A ← {qs‿gs𝕊v:
     4│      1⊸=◶{𝕊𝕩: ui ← 0 L gs
 22430│        v M˜ {𝕊⟨⟩:ui; (M˜´M⟜ui⟜M˜M´) L⟜g.swap¨ 𝕩} T qs (⌽∘⊢∾¬∘∊/⊣)˜ ↕qb
   366│      }‿(v M˜ gs L˜ ⊑qs) ≠qs}
     5│    »⊸<∨` 0>r-`>+○(×˜)˝ wf A´ ⌽sc
      │  }

Try running the simulation in the BQN repl and explore it! If you are an Emacs user, the org-mode computational notebook in the blog's repository provides the best experience.

Footnotes:

Although conceptually straightforward, a Hilbert space of size 2⋆n makes this type of simulation a true computational challenge. For an efficient implementation, see arXiv:1710.05867.

This post's title is a playful recursive acronym that employs quantum computing terminology, without any specific significance beyond that.

I optimized it up to this number, but I wasn't targeting the Kolmogorov complexity.

⁴

Hilbert's radio address in 1930.

⊑∘∞