Towards a Linear-Algebraic Hypervisor
arXiv:2604.12902v1 [cs.PL] 14 Apr 2026
Breandan Considine Abstract
2
Many techniques in program synthesis, superoptimization, and array programming require parallel rollouts of generalpurpose programs. GPUs, while capable targets for domainspecific parallelism, are traditionally underutilized by such workloads. Motivated by this opportunity, we introduce a pleasingly parallel virtual machine and benchmark its performance by evaluating millions of concurrent array programs, observing speedups up to 147× relative to serial evaluation.
We will consider a parallel RASP variant loosely equivalent to the PRAM [5] model, except each machine stores its program in memory and all machines occupy constant space. Concretely, we lift the transition function to a vectorized transition map, Φ : 𝐶 𝑑 → 𝐶 𝑑 whose iterates are given by Φ𝑖+1 = Φ ◦ Φ𝑖 , so the resulting configuration after 𝑖 steps is c𝑖 = Φ𝑖 (c0 ). Thus, for a vector of programs P and inputs x, we will write c𝑖 (P, x) := Φ𝑖 c0 (P, x) for the evolution of 𝑛 independent machines, each with static memory. Since memory is statically allocated, each stepwise transition is a Boolean function on finitely many words, hence Φ admits a unique multilinear polynomial representation over F2 [8].
1
Introduction
An abstract machine is a sequential model of computation for analyzing the complexity of algorithms, which, unlike the Turing machine, more closely characterizes the operational behavior of modern computers. Elgot and Robinson [4] give an early example of such a machine, known as the randomaccess stored program (RASP) machine, equipped with a small set of atomic instructions and a shared address space storing both program and data. Cook and Reckhow [3] later present a simplified RASP and introduce an equivalent model which isolates its program from its data in memory, known as the random-access machine (RAM). An implementation of an abstract machine is known as a virtual machine (VM) and a hypervisor is a virtualization environment that simulates multiple VMs on a single physical machine, such as a GPU. 1.1
RASP model
A RASP machine is a 4-tuple ℛ = ⟨𝐶, 𝑐 0, →, 𝐹 ⟩ consisting of a configuration space, 𝐶 = N×Z×ZN ×Z∗ ×Z∗ , which we will write as ⟨𝑖, 𝑎, 𝑀, 𝑢, 𝑦⟩ : 𝐶, with 𝑖 being the instruction counter, 𝑎 being the accumulator, 𝑀 being the memory, 𝑢 being the unread input, and 𝑦 being the output. A RASP program is a vector 𝑃 : (Z × Z)𝑚 that is packed with an input 𝑥 : Z∗ , into ′ an initial configuration, 𝑐 0 (𝑃, 𝑥) as ⟨0, 0, 𝑀0..2𝑚 ← 𝑃, 𝑥0, 𝜀⟩. Letting ⟨𝑜, 𝑗⟩ = ⟨𝑀𝑖 , 𝑀𝑖+1 ⟩ denote the opcode and operand respectively, the one-step transition function → : 𝐶𝐶 is ⟨𝑖 + 2, 𝑗, 𝑀, 𝑢, 𝑦⟩, ⟨𝑖 + 2, 𝑎 + 𝑀 𝑗 , 𝑀, 𝑢, 𝑦⟩, ⟨𝑖 + 2, 𝑎 − 𝑀 𝑗 , 𝑀, 𝑢, 𝑦⟩, ⟨𝑖 + 2, 𝑎, 𝑀 𝑗 ← 𝑎, 𝑢, 𝑦⟩, ⟨𝑖, 𝑎, 𝑀, 𝑢, 𝑦⟩ → ⟨𝑗, 𝑎, 𝑀, 𝑢, 𝑦⟩, ⟨𝑖 + 2, 𝑎, 𝑀, 𝑢, 𝑦⟩, ⟨𝑖 + 2, 𝑎, 𝑀 𝑗 ← 𝜎, 𝑣, 𝑦⟩, ⟨𝑖 + 2, 𝑎, 𝑀, 𝑢, 𝑦 · 𝑀 𝑗 ⟩, ⟨𝑖, 𝑎, 𝑀, 𝑢, 𝑦⟩,
𝑜 = 1LOD, 𝑜 = 2ADD, 𝑜 = 3SUB, 𝑜 = 4STO, 𝑜 = 5BPA ∧ 𝑎 > 0, 𝑜 = 5BPA ∧ 𝑎 ≤ 0, 𝑜 = 6RD ∧ 𝑢 = 𝜎𝑣, 𝑜 = 7PRI, otherwise.
The RASP machine’s halting configurations, 𝐹 ⊆ 𝐶, are characterized by the fixed points, 𝐹 = {𝑐 ∈ 𝐶 | 𝑐 → 𝑐}.
2.1
Method
Word RASP
Fix 𝑤, 𝑚, 𝑛, ℓ, 𝑠 ∈ N, and assume a fixed-length word, W = F2𝑤 . Then, for a program 𝑃 : W2𝑚≤𝑛 , define the word RASP as ℛ (𝑤,𝑛,ℓ,𝑠 ) = ⟨𝐶 𝑤 , 𝑐 0, →, 𝐹 ⟩ with a finite configuration space, 𝐶 𝑤 = W × W × W𝑛 × Wℓ+1 × W𝑠+1 = F2𝑤 (𝑛+ℓ+𝑠+4) , 𝑐 = ⟨𝑖 , 𝑎 , 𝑀 , 𝑢 , 𝑦⟩ where the initial configuration on input 𝑥 is given by 𝑐 0 (𝑃, 𝑥) = 0𝑤 , 0𝑤 , 𝑃 0𝑤𝑛−2𝑚 , 𝑥0ℓ𝑤− |𝑥 | , 0𝑠+1 . 𝑤 We write ⊕, ⊗ : W × W → W for addition and multiplication (mod 2𝑤 ) and assume a circular memory addressing scheme, ⟨𝑜, 𝑗⟩ := ⟨𝑀 𝑖 mod 𝑛 , 𝑀 (𝑖 ⊕1) mod 𝑛 ⟩. Now, let 𝛿 𝑗𝑘 be the Kronecker delta and define the operators, 𝑛−1 𝑀 𝑗 ← 𝑎 := 𝛿 𝑗𝑘 𝑎 + (1 − 𝛿 𝑗𝑘 )𝑀𝑘 𝑘=0 , 𝑢 • := 𝑢 min(𝑢0 +1,ℓ ) , 𝑢 ⊲ := 𝑢 0 ← min(𝑢 0 + 1, ℓ), ( (𝑦 𝑦0 +1 ← 𝑧)0 ← 𝑦0 + 1, 𝑦0 < 𝑠, 𝑦 :: 𝑧 := 𝑦, otherwise. The one-step transition operator takes ⟨𝑖, 𝑎, 𝑀, 𝑢, 𝑦⟩ to ⟨𝑖 ⊕ 2𝑤 , 𝑗, 𝑀, 𝑢, 𝑦⟩, 𝑜 = 1LOD, ⟨𝑖 ⊕ 2𝑤 , 𝑎 ⊕ 𝑀 𝑗 mod 𝑛 , 𝑀, 𝑢, 𝑦⟩, 𝑜 = 2ADD, ⟨𝑖 ⊕ 2 , 𝑎 ⊗ 𝑀 , 𝑀, 𝑢, 𝑦⟩, 𝑜 = 3MUL, 𝑤 𝑗 mod 𝑛 ⟨𝑖 ⊕ 2 , 𝑎, 𝑀 ← 𝑎, 𝑢, 𝑦⟩, 𝑜 = 4STO, 𝑤 𝑗 mod 𝑛 ⟨𝑗, 𝑎, 𝑀, 𝑢, 𝑦⟩, 𝑜 = 5BNZ ∧ 𝑎 ≠ 0, ⟨𝑖 ⊕ 2𝑤 , 𝑎, 𝑀, 𝑢, 𝑦⟩, 𝑜 = 5BNZ ∧ 𝑎 = 0, ⟨𝑖 ⊕ 2𝑤 , 𝑎, 𝑀 𝑗 mod 𝑛 ← 𝑢 •, 𝑢 ⊲ , 𝑦⟩, 𝑜 = 6RD ∧ 𝑢 0 < ℓ, 𝑖 ⊕ 2𝑤 , 𝑎, 𝑀, 𝑢, 𝑦 :: 𝑀 𝑗 mod 𝑛 , 𝑜 = 7PRI, ⟨𝑖, 𝑎, 𝑀, 𝑢, 𝑦⟩, otherwise. Once again, the word RASP machine’s halting configurations are characterized by the fixed points, 𝐹 = {𝑐 ∈ 𝐶 𝑤 | 𝑐 → 𝑐}.
Breandan Considine
Let Φ𝑤 : 𝐶 𝑤 → 𝐶 𝑤 denote the resulting coordinatewise polynomial, so 𝑐 → 𝑐 ′ ⇐⇒ Φ(𝑐) = 𝑐 ′ . Alternatively, this can be written as a total deterministic transition function, Ê Ê Φ𝑤 (𝑐) = 𝜂𝜓 (𝑐) Φ𝜓 (𝑐), 𝜂𝜓 (𝑐) = 1, 𝜓
𝜓
where 𝜓 ranges over the guarded instruction cases, 𝜂𝜓 is a multilinear indicator, and Φ𝜓 is the associated branch update. 2.2
Heapless array lowering
3
Evaluation
We implement VMs for the word RASP on (1) Kotlin/JVM and (3) Nvidia CUDA with parameters 𝑤 = 32, 𝑛 = 250, ℓ = 10, 𝑠 = 2, 𝜏max = 106 . Then we evaluate each implementation on 𝑑 = 8 × 106 random array programs sampled uniformly from the space of valid code snippets of exactly length 𝐿 = 100 using Considine’s word sampler [2]. In Fig. 1a, we measure runtime across the JVM (running on an Apple M4 Max) and CUDA platforms (on the Nvidia A10G and B200 GPUs), reporting wall-clock time to evaluate up to 𝜏max VM steps.
Consider a small array language with the following terms: (b) Halting distribution
(a) Average runtime
F ::= fun f0 ( ipt : W ^ N ) -> W ^ N { B }
106
N ::= 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 Time (s)
B ::= S | S ; B | hlt S ::= P = E | ife G { B } { B } | whl G { B } Z ::= 0 | N G ::= ipt [ Z ] | scr [ Z ] | opt [ Z ]
M4 (×1)
M4 (×16)
B200 (×1)
B200 (×8)
102 101
102
𝜏ℎ ∈ [102 , 106 − 1] →
1 2 3 4 5 6 7 8
10 20 30 40 50 60 70 80 9 10 0 0+
E ::= G + G | G * G | G + Z | G * Z | G | Z
104
0
0
P ::= scr [ Z ] | opt [ Z ]
Reserve disjoint memory regions 𝑝 0...𝑙 −1 , 𝑞 0...𝑒 −1 , and 𝑟 0...𝛾 −1 for arrays ipt : W𝑙 ≤ℓ , opt : W𝑒 ≤𝑠 and scr : W𝛾 ≤𝜇 , with (5 + ℓ + 𝑠 + 𝜇 + 2𝑚) ≤ 𝑛 and interpret array indexing modulo declared length, so J𝛼K𝑧 := 𝑀2𝑚+J𝛼 K+(𝑧 mod |𝛼 | ) . We will now assign a small-step operational semantics to this language:
106 ≤ 𝜏ℎ ≈ ∞ →
A10G
Count
103
Total VMs (×106 )
Halting steps (𝜏ℎ )
Figure 1. Wallclock and simulation runtime of 8m programs. In Fig. 1b, we sample 𝑑 = 8 × 106 programs of length 𝐿 = 100 and estimate the halting probability before 𝜏max = 106 steps, 𝑤·arity(𝑃 ′ ) Pr Φ𝜏𝑤max 𝑐 0 (𝑃 ′ ⇒ 𝑃, 𝑥 ∼ F2 ∈ 𝐹 𝐿 = |𝑃 ′ | . ′ 𝑃 ∼𝒫𝐿
Γ ⊢𝐵 ⇒𝐶 −1 HLT ⊢ fun f0 ( ipt : I ^ N ) -> I ^ 𝑙 { 𝐵 } ⇒ 𝐶 PRI 𝑞𝑖 𝑙𝑖=0
Γ ⊢ 𝛼 [𝑧 ] ⇒𝑎 LOD 0 ADD J𝛼 K𝑧 Γ ⊢ 𝐼 1 ⇒𝑎 𝐽1
Get
Γ ⊢ 𝐼 2 ⇒𝑎 𝐽2
Γ ⊢ 𝐼 1 + 𝐼 2 ⇒𝑎 𝐽1 STO 𝜈 𝐽2 ADD 𝜈 Γ ⊢ 𝐼 ⇒𝑎 𝐽
⊕
Γ ⊢ 𝐸 ⇒𝑎 𝐴 Γ ⊢ 𝛼 [𝑧 ] = 𝐸 ⇒ 𝐴 STO J𝛼 K𝑧 Γ ⊢ 𝐼 1 ⇒𝑎 𝐽1
Γ ⊢ 𝐵1 ⇒ 𝐶1
Put
4
Γ ⊢ 𝐼 2 ⇒𝑎 𝐽2
Γ ⊢ 𝐵2 ⇒ 𝐶2
Γ ⊢𝐵 ⇒𝐶
Γ ⊢ whl 𝐼 {𝐵 } ⇒ ℓℎ : 𝐽 BNZ ℓ𝑏 JMP ℓ𝑒 ℓ𝑏 : 𝐶 JMP ℓℎ ℓ𝑒 : Γ ⊢ 𝐵 ⇒ 𝐶𝐵 Γ ⊢ 𝑆 ⇒ 𝐶𝑆 Seq Γ ⊢ 𝐵 ; 𝑆 ⇒ 𝐶 𝐵 𝐶𝑆
We plot the distribution of halting programs by their halting number 𝜏ℎ , i.e., the least 𝜏 such that Φ𝜏𝑤 (𝑐 0 ) = Φ𝜏+1 𝑤 (𝑐 0 ).
Γ ⊢ 𝐼 1 × 𝐼 2 ⇒𝑎 𝐽1 STO 𝜈 𝐽2 MUL 𝜈
Γ ⊢ ife 𝐼 {𝐵 1 } {𝐵 2 } ⇒ 𝐽 BNZ ℓ𝑡 𝐶 2 JMP ℓ𝑒 ℓ𝑡 : 𝐶 1 ℓ𝑒 : Γ ⊢ 𝐼 ⇒𝑎 𝐽
Fun
⊗
Ife
While
−1 HLT Γ ⊢ hlt ⇒ PRI 𝑞𝑖 𝑙𝑖=0
Halt
We write Γ ⊢ 𝐸 ⇒ 𝐴 when the instruction sequence 𝐴 implements 𝐸, and Γ ⊢ 𝐸 ⇒𝑎 𝐴 when 𝐴 implements 𝐸 and leaves value 𝑎 in the accumulator. Let 𝜈 be a single reserved memory address and finally, desugar JMP ℓ as LOD 1 BNZ ℓ , with HLT being any reserved invalid instruction pair whose execution would trigger a fixed-point halt (e.g., 0 0 ). 2.3
Implementation
VMs are implemented using a mainly branchless strategy to reduce warp divergence on the GPU, and host threads are dispatched to VMs using a simple round-robin scheduler. A short supplemental, including hypervisor pseudocode and platform-specific details, may be found in Appendix A.
Related work
Prior work has explored compiling to transformers [10], though not intended as a practical compiler per se. The most relevant systems literature is Vectorvisor [6], which targets WebAssembly IR, requiring a much larger set of opcodes. Exhaustive search has played a key role in recent empirical mathematics projects [1, 7] which use distributed computing to filter for constructions with favorable properties, but do not explicitly leverage SIMD architectures. We specifically target a much smaller IR and heap space, enabling us to run millions of concurrent VMs. The RASP model’s simplicity, low memory footprint, and proximity to realistic assembly makes it suitable for algorithmic analysis, program synthesis, empirical research, and a variety of pedagogical applications.
5
Conclusion
We have presented a multilinear RASP VM, illustrated its viability by lowering a heapless array language, and used it to estimate the probability of sampling halting programs, demonstrating significant speedups over naïve evaluation. In future work, we intend to use it to accelerate the search for straightline programs for fast matrix multiplication and context-free grammar parsing. Source code and data for all experiments will be released here: https://github.com/breandan/lavm
Towards a Linear-Algebraic Hypervisor
References
B
[1] Justin Blanchard, Daniel Briggs, Konrad Deka, Nathan Fenner, Yannick Forster, Georgi Georgiev, Matthew L House, Rachel Hunter, Maja Kądziołka, Pavel Kropitz, et al. 2025. Determination of the fifth Busy Beaver value. arXiv preprint arXiv:2509.12337 (2025). [2] Breandan Considine. 2025. A word sampler for well-typed functions. arXiv:2512.01036 [cs.PL] https://arxiv.org/pdf/2512.01036.pdf [3] Stephen A Cook and Robert A Reckhow. 1972. Time-bounded random access machines. In Proceedings of the fourth annual ACM symposium on Theory of computing. 73–80. [4] Calvin C Elgot and Abraham Robinson. 1964. Random-access storedprogram machines, an approach to programming languages. Journal of the ACM (JACM) 11, 4 (1964), 365–399. [5] Steven Fortune and James Wyllie. 1978. Parallelism in random access machines. In Proceedings of the tenth annual ACM symposium on Theory of computing. 114–118. [6] Samuel Ginzburg, Mohammad Shahrad, and Michael J Freedman. 2023. VectorVisor: a binary translation scheme for throughput-oriented GPU acceleration. In 2023 USENIX Annual Technical Conference. 1017–1037. [7] Manuel Kauers and Jakob Moosbauer. 2023. Flip graphs for matrix multiplication. In Proceedings of the 2023 International Symposium on Symbolic and Algebraic Computation. 381–388. [8] Elchanan Mossel, Ryan O’Donnell, and Rocco P Servedio. 2003. Learning juntas. In Proceedings of the thirty-fifth annual ACM symposium on Theory of computing. 206–212. [9] NVIDIA Corporation. 2026. CUDA Programming Guide. https://docs. nvidia.com/cuda/cuda-programming-guide/index.html Release 13.2. [10] Gail Weiss, Yoav Goldberg, and Eran Yahav. 2021. Thinking like transformers. In International Conference on Machine Learning. PMLR, 11080–11090.
Listed below are a few of the busiest beavers discovered after conducting a ∼ 20 minute search with the JVM running on an M4 (×16) architecture using the grammar from Sec. 2.2, with the ipt, opt and scr parameters all initialized to 0s.
A
Implementation details
We launch 𝑊 persistent workers, where 𝑊 is the maximum number of resident threads, and stripe the global configuration vector c, so that worker 𝑔 visits each VM in its stripe in round-robin order. Each visit advances one live VM for a bounded epoch 𝑞, yielding a static schedule that realizes the lifted dynamics up to a finite time horizon, Φ𝜏max : 𝐶 𝑑𝑤 → 𝐶 𝑑𝑤 . Algorithm 1 Striped round-robin RASP hypervisor 1: inputs: Global config c0 ∈ 𝐶 𝑑 , epoch 𝑞, step budget 𝜏max 2: c ← c0 , 𝑊 ← MaxResidentThreads() 3: for 𝑔 = 0 to 𝑊 − 1 in parallel do ⊲ Persistent workers 4: for 𝑟 = 0 to ⌈𝜏max /𝑞⌉ − 1 do ⊲ Total rounds 5: for 𝑘 = 0 to ⌈𝑑/𝑊 ⌉ − 1 do ⊲ Walk worker stripe 6: 𝑗 ← 𝑔 + 𝑘𝑊 ⊲ Current VM index 7: if 𝑗 < 𝑑 and c[ 𝑗] ∉ 𝐹 then ⊲ Skip halted VMs 8: 𝑅 ← LoadRegs(c[ 𝑗]) 9: for 𝑡 = 1 to 𝑞 do ⊲ Bounded epoch 10: if 𝑅 ∈ 𝐹 then break 11: 𝑅 ← Φ𝑤 (𝑅) ⊲ One-step transition 12:
c[ 𝑗] ← StoreRegs(𝑅)
13: return c
Since CUDA is limited to at most 255 registers per thread [9], we constrain VM memory to 𝑛 ≲ 250 32-bit words, however this constraint can be relaxed, VRAM permitting, to 1 kB or higher, albeit at the cost of increased memory traffic. Assuming 256 resident threads per SM, an Nvidia B200 can sustain about 37,888 resident threads at peak theoretical occupancy.
Heapless busy beavers (𝐿 = 120)
BB #1 (𝜏ℎ = 1549) fun f0(ipt: W ^ 4) -> W ^ 1 { opt[0] = 4; whl opt[0] { opt[0] = opt[0] + opt[0]; scr[3] = opt[0]; whl opt[0] { whl opt[0] { scr[7] = scr[2] * 7; opt[0] = 0; scr[1] = 9 } }; opt[0] = scr[2] + 0; opt[0] = scr[3] + scr[5] } }
BB #2 (𝜏ℎ = 1272) fun f0(ipt: W ^ 2) -> W ^ 1 { opt[0] = opt[0] + 6; whl opt[0] { scr[3] = opt[0] + 0; ife ipt[0] { hlt } { scr[8] = 6; opt[0] = ipt[0] + 3; scr[0] = ipt[1] * 2; opt[0] = scr[3] * 6; scr[6] = scr[0] * 3; opt[0] = opt[0] } } }
BB #3 (𝜏ℎ = 1255) fun f0(ipt: W ^ 3) -> W ^ 1 { opt[0] = 1; whl opt[0] { opt[0] = opt[0] + ipt[0]; ife ipt[2] { hlt } { scr[4] = 6; opt[0] = opt[0] + opt[0] }; scr[6] = scr[8] * 2; whl ipt[1] { scr[4] = 3; opt[0] = 8 }; scr[0] = opt[0] + 9 } }