Jianming Tong∗
Jingtian Dang∗
Simon Langowski
[email protected] Georgia Institute of Technology Atlanta, Georgia, USA
[email protected] Georgia Institute of Technology Atlanta, Georgia, USA
[email protected] Massachusetts Institute of Technology Cambridge, Massachusetts, USA
Tianhao Huang
Asra Ali
Jeremy Kun
[email protected] Massachusetts Institute of Technology Cambridge, Massachusetts, USA
[email protected] Google Austin, Texas, USA
[email protected] Google Portland, Oregon, USA
Jevin Jiang
Srinivas Devadas
Tushar Krishna
[email protected] Google Sunnyvale, California, USA
[email protected] Massachusetts Institute of Technology Cambridge, Massachusetts, USA
[email protected] Georgia Institute of Technology Atlanta, Georgia, USA
Multi-Scalar Multiplication Number Theory Transformation Elliptic Curve Multiplication Modular Multiplication over Large Prime Finite Field 256 bits Scalar +ω +ω2 x y Point -ω2 + zz +ω2 zzz 377 bits NTT MSM 377 bits
Abstract Zero-knowledge proof (ZKP) provers remain costly because multiscalar multiplication (MSM) and number-theoretic transforms (NTTs) dominate runtime as they need significant computation. AI ASICs such as TPUs provide massive matrix throughput and SotA energy efficiency. We present MORPH, the first framework that reformulates ZKP kernels to match AI-ASIC execution. We introduce Big-𝑇 complexity, a hardware-aware complexity model that exposes heterogeneous bottlenecks and layout-transformation costs ignored by Big-𝑂. Guided by this analysis, (1) at arithmetic level, MORPH develops an MXU-centric extended-RNS lazy reduction that converts high-precision modular arithmetic into dense lowprecision GEMMs, eliminating all carry chains, and (2) at dataflow level, MORPH constructs a unified-sharding layout-stationary TPU Pippenger MSM and optimized 3/5-step NTT that avoid on-TPU shuffles to minimize costly memory reorganization. Implemented in JAX, MORPH enables TPUv6e8 to achieve up-to 10× higher throughput on NTT and comparable throughput on MSM than GZKP. Our code: https://github.com/EfficientPPML/MORPH.
1
MORPH
arXiv:2604.17808v1 [cs.AR] 20 Apr 2026
Enabling AI ASICs for Zero Knowledge Proof
Dataflow Optimizations: Maximize Parallelism and Minimize Layout Transformation Arithmetic Optimizations: MXU-Centric RNS Lazy Reduction for Prime Field Mul. XLA Compiler AI ASICs E.g.
Google TPU
AWS Trainium
Figure 1: MORPH’s deployment flow with dataflow and arithmetic optimizations to accelerate ZKP on TPU. Both MSM and NTT expose abundant parallelism at practical ZKP sizes, making parallel hardware promising for acceleration. Among these platforms, AI ASICs (Google TPU [21], AWS Trainium [4] etc.) offer extreme compute density and energy efficiency, significantly outperforming general-purpose accelerators at scale [35]. Architecturally, a TPU contains giant heterogeneous cores, including a matrix multiplication engine (MXU) and a vectorized processing engine (VPU). Each MXU/VPU contains orders of magnitude (1024/16×) more multiply-accumulates (MACs) than one corresponding core in a GPU [8]. All MACs in a MXU/VPU execute instruction in the lock-step (SIMD) to substantially reduce per-operation control overhead and improve energy efficiency. This paper studies how to systematically repurpose AI ASICs, particularly Google’s TPUs [21], for ZKP kernels (MSM/NTT) and achieve higher energy efficiency and SoTA throughput at scale. Existing ZKP acceleration spans GPUs [14, 20, 27, 28], FPGAs [7, 26, 29–32] and ASICs [33, 38, 42]. While these systems hand-craft MSM/NTT algorithms for their respective platforms, they provide no principled way to quantify how far an MSM/NTT algorithm is from the platform’s performance limits with detailed reasoning. Some works provide Big-O complexity [22], yet we find a surprising phenomenon: the fastest GPU/TPU algorithms intentionally use the algorithm with higher Big-O to expose massive parallelism and thus achieve higher throughput [35]. Big-O alone is insufficient for modern heterogeneous parallel hardware. Furthermore, Big-O cannot capture layout transformation cost, which becomes the dominant bottleneck on AI ASICs. When we
Introduction
Zero-knowledge proofs (ZKPs) have evolved from a theoretical construct into a core building block for practical verifiable computation and scalable blockchains. They allow services to outsource computation while keeping data private and verification cheap, but the prover remains prohibitively expensive, e.g., generating a proof for ImageNet ViT requires nearly one hour [41]. Across widely deployed protocols [17, 18], prover runtime is dominated by multi-scalar multiplication (MSM) and the number-theoretic transform (NTT), which commonly account for ≈70% and 20–30% of total latency, respectively [20, 42]. ∗ Both authors contributed equally to this research.
This work is licensed under a Creative Commons Attribution 4.0 International License. DAC ’26, Long Beach, CA, USA © 2026 Copyright held by the owner/author(s). ACM ISBN 979-8-4007-2254-7/2026/07 https://doi.org/10.1145/3770743.3804402 1
• BIG-𝑇 Complexity: A hardware-aware complexity metric that captures heterogeneous compute spans and layout transformation, enabling principled analysis of MSM/NTT designs on AI 1 MORPH: Matrix-Oriented Reformulation of ZKP for Parallel Heterogeneous AI ASICs
2
10~200 GB/s PCI-E
HBM 10~100 GB
CPU 1~10 TOPS
~1000 GB/s VMEM ~0.1 GB TPU
port the SotA MSM and NTT algorithms for GPU/FPGA/ZKP-ASIC to TPUs with similar overall throughput, they slow down by 30× in our evaluation. This is because a TPU only operates on coarsergrained contiguous data (a 4KB vector register, VReg), and fine-grained data gathering in SotA MSM/NTT algorithms requires expensive shuffles and transposes to move data into contiguous VRegs. These layout-induced stalls are invisible under the traditional Big-O model. To address this gap, we introduce Big-𝑇 notation, a platform-aware complexity model that measures the sequential bottleneck across heterogeneous pipelined compute units and the latency of layout transformation, providing an asymptotic lower bound for parallel execution on AI ASICs. Big-𝑇 exposes two structural inefficiencies in SotA MSM/NTT algorithms on TPUs. Arithmetic-level Challenge: ZKP requires modular arithmetic over large prime fields (e.g., 256/377/753 bits), but hardware natively supports much lower precision (e.g., 8/32-bit). Because prime fields lack a natural RNS factorization, the SotA algorithm [6, 27] performs high-precision arithmetic via digit decomposition and a chain of carry propagation, achieving <1% compute utilization on TPUs. Dataflow-level Challenge: The SotA dataflows incur prohibitive communication and storage overheads. In presorted Pippenger MSM [20, 27], sorting is distributed across many GPU CUDA cores. Directly porting this strategy to TPU introduces prohibitive interTPU communication. Layout-invariant 3-step NTT [35] requires parameters that scale linearly with problem size, resulting in outof-memory (OOM) failures at large scale. In resolving these inefficiencies, we propose MORPH1 , the first framework enabling TPU as ZKP accelerator with arithmetic and dataflow level optimizations, as shown in Fig. 1. Arithmetic-level Solution: Since a prime field cannot be represented in a Residue Number System (RNS), MORPH adopts [19, 24] to encapsulate modular multiplication in a prime field in a larger non-prime finite field that has a RNS representation. This enables high-precision multiplication to be expressed as independent 32-bit vector multiplications with no carry propagation. Then, MORPH reforms the modular reduction (RNS reconstruction → modular reduction → RNS decomposition) as low-precision matrix multiplication and parallel vectorized operations to be accelerated by TPU’s MXU and VPU, achieving up-to 90× speedup. Dataflow-level Solution: To tackle communication and layout reorganization, MORPH proposed LS-PPG to schedule communicationfree workload dimensions across devices to avoid inter-TPU communication, ensures layout invariant within individual TPU. To mitigate parameters storage overflow, MORPH’s 5-step NTT recursively partitions workloads to reduce parameter size while preserving the same computational complexity. Together, MORPH enables TPUv6e-8 to achieve up-to 10× NTT throughput and 1.2× MSM throughput than GZKP [27] (NVIDIA V100). MORPH makes TPU the SotA throughput machine for high-precision NTT, a promising platform for ZKP acceleration. Our contributions are:
32 VReg 4 KB 8 Compute
Max Parallelism Datatype Requirement
Lane 0 0.125~1 MB
Lane 1
Lane 127
~10 TB/s 32-bit regs VPU SIMD 1~10 TOPS
Shuffle 𝑃𝐴𝑅S = 4096 u32 (8, 128)
"SIMD" Coarse-Grained Control
128 MXU Systolic Array (MatMul) 100~1000 TOPS
Transpose 𝑃𝐴𝑅T = 4096 u32 (8, 128)
XLU Layout Transform
Scalar Unit SReg, SMEM For control
MXU 𝑃𝐴𝑅MXU = 131072 bf16/i8 (8×128×128 MatMul)
VPU 𝑃𝐴𝑅VPU = 2048 u32 2048
Figure 2: TPU programming model. All compute units operate on VRegs with (8 × 128) 32-bit registers each. Performance hinges on data being laid out such that tiles needed for computation can be loaded directly from VRegs. Values from [21]. ASICs and guiding the construction of layout-stationary pippenger for MSM and 5-step NTT, giving 5.8× and 1.6× speedup. • MXU-centric RNS Lazy Reduction: A matrix-oriented modular reduction for big integer multiplication in prime fields, which replaces sequential carry-propagation arithmetic with MXUaccelerated RNS-lazy reduction, achieving up-to 90× speedup over the SotA radix decomposed Montgomery reduction. • LS-PPG and 5-step Layout Invariant NTT: Both avoid interTPU communication and minimize intra-TPU layout reorganization, delivering up-to 5.8× speedup over the SoTA MSM algorithm[20] on TPU, and reduce parameters storage at scale.
2 Background and Related Work 2.1 Tensor Processing Unit (TPU) The TPU [21] is Google’s AI ASIC designed to accelerate machine learning training and inference with high energy efficiency. All on-chip computation is organized around the Vector Register (VReg) abstraction. Each VReg consists of 8×128 of 32-bit values executed in lock-step SIMD. All compute units on a TPU, including the VPU, MXU, and XLU, consume and produce data strictly at VReg granularity, as shown in Fig. 2. This uniform but coarse layout enables high efficiency but restricts flexibility: workloads whose shapes do not align with the fixed (8, 128) structure incur layout transformation because only full VRegs can be processed each cycle. • Vector Processing Unit (VPU). The VPU performs 32-bit SIMD arithmetic on VRegs. E.g., TPUv4 contains a dual-issue ALU to operate two VRegs concurrently, yielding a peak parallelism of 𝑃VPU = 2048 32-bit operations per cycle. Vectors with length not a multiple of (8, 128) cannot fully utilize VPU’s available parallelism. • Matrix Multiplication Unit (MXU). The MXU in TPUv4 is a large systolic array (128×128) for int8 and bf16 matrix multiplication. It preloads weight tiles from VRegs and streams activation tiles from VRegs. When normalized to the same datatype, the MXU achieves roughly 16× higher peak throughput than the VPU. Full efficiency requires workloads to conform to the MXU’s minimum effective tile of 8×128×128. Larger gives more weights reuse while smaller or non-aligned matrix shapes leave MXU under-utilized. • Layout Transformation Unit (XLU). XLU performs shuffles, transposes, broadcasts, and reductions of VRegs. XLU is efficient for transformations at VRegs granularity, but costly for element-wise
is split into 𝐾 = 𝑆 𝐵𝐿 /𝑐 windows of 𝑐 bits. We adopt presorting PPG as our dataflow baseline, termed as “Presort-PPG”. Butterfly NTT and layout invariant 3-step NTT. SotA NTT acceleration on CPUs [11, 23], GPUs [14, 27], FPGAs [32] and ASICs [? ] implements recursive Cooley-Tukey NTT with minimal computation complexity 𝑂 (𝑁𝑙𝑜𝑔𝑁 ), commonly called as “butterfly NTT" in Fig. 5a. However, each stage performs fine-grained shuffles that are smaller than vector-register width, leading to costly layout transformations, making TPUs inefficient. Recent TPU work [35] proposes a layout-invariant 3-step NTT achieving 𝑂 (𝑁 3/2 ) arithmetic with zero layout transformation, breaking the NTT throughput record of the butterfly NTT on GPUs [12, 16] on lower-precision 𝑀 < 32 bits and lower degree (𝑁 ≤ 217 ). It reforms an 𝑁 -input NTT into an (𝑅, 𝐶) matrix (𝑁 = 𝑅 × 𝐶) such that an NTT is reformed into matrix and vectorized multiplication for TPU acceleration, as shown in Fig. 5b. We adopt 3-step NTT as dataflow baseline.
permutations. Shuffling 𝑁 individual elements may cost 𝑁 cycles. • HBM and the host CPU serve as the off-chip data sources, but their bandwidth is typically an order of magnitude lower than on-chip VReg bandwidth. Consequently, high-performance TPU kernels favor on-chip data reuse and layout transformation.
2.2
Performance Metric and Analysis
TPU’s heterogeneous compute units, each with different parallelism and tiling constraints, make classical models such as Big 𝑂 and the work–span model inadequate. These models measure total computational work or critical-path latency under the assumption of a homogeneous sequential machine like a CPU, and therefore cannot capture VReg granularity, systolic-array tiling, or the cost of on-chip layout transformations. Further, Roofline analysis[39] and trace-based profiling quantify how far execution is from hardware limits, but they still do not explain why performance is low.
2.3
3
Zero-knowledge Proof (ZKP) Acceleration
Method
This section first defines the Big 𝑇 notation, then uses it to analyze arithmetic and dataflow inefficiencies of the SotA MSM/NTT algorithm on TPUs, finally showing how MORPH resolves them.
2.3.1 Background on Kernel and Arithmetic in ZKP. Í𝑁 −1 𝑆𝑛 ⊛ Multi-Scalar Multiplication. MSM [18] computes 𝑛=0 𝑃𝑛 where each scalar 𝑆𝑛 is 𝑆 𝐵𝐿 bits and each point 𝑃𝑛 lies on an elliptic curve that supports point addition (PADD, ⊕). A scalar–point multiplication (⊛) is realized as repeated PADDs for 𝑆𝑛 times. Number Theory Transformation. NTT [18] converts a polynomial of degree 𝑁 into its frequency-domain representation by evaluating it at 𝑁 distinct roots of unity in a finite field. Given an input Í𝑁 −1 vector x of degree 𝑁 , NTT computes X𝑘 = 𝑛=0 𝑥𝑛 × 𝜔 𝑘𝑛 𝑁 mod 𝑀, 𝑘 ∈ [0, 𝑁 ), where 𝜔 𝑁 is a primitive 𝑁 -th root of unity. Large prime field. Modern ZKP systems [18] rely on arithmetic over a finite field F𝑀 with a large prime modulus 𝑀 (256/377/753 bits). In both elliptic-curve MSM and polynomial-commitment NTT, all additions, multiplications are executed modulo large prime 𝑀. 2.3.2 Related Work: Arithmetic Implementations. ZKPs operate over large prime fields F𝑀 with log2 (𝑀) ≥ 256, exceeding native 32-bit word size on GPUs and TPUs. SotA GPU implementations [20, 27] therefore represent each field element as a vector of 𝐷 32-bit “digits”. A single field multiplication requires two stages: (1) High-precision multiplication, which performs 𝑂 (𝐷 2 ) 32-bit multiplications followed by a 𝑂 (𝐷) sequential carry propagation of length 𝐷; and (2) Montgomery reduction, which reduces the 2𝐷-digit intermediate product back to 𝐷 digits, incurring another 𝑂 (𝐷 2 ) 32-bit multiplications. We use this SotA radix Montgomery reduction as our arithmetic baseline, termed as radix Mont. RNS represents an integer 𝑥 by residues (𝑥 mod 𝑞𝑖 ) over coÎ prime moduli {𝑞𝑖 }, where modulus 𝑄 = 𝑖 𝑞𝑖 . RNS could lower 𝑂 (𝐷 2 ) down to 𝑂 (𝐷) when performing modular multiplication with respect to 𝑄. But the prime modulus 𝑀 cannot be factored into multiple residues, making such reduction inapplicable to ZKP. 2.3.3 Related Work: Dataflow Optimizations. Presorting Pippenger (PPG) - MSM In prior work of MSM acceleration on GPU [20, 28], FPGA [7, 26, 29–31] and ASICs [13, 25, 38, 40, 42], Presorting PPG is the SotA MSM algorithm. Its key idea is to reduce the total number of point additions by first grouping points whose scalars share the same value, summing those points once, and then multiplying the accumulated result by that scalar value. To increase the likelihood of such sharing (“collisions”), each scalar
3.1
Big-𝑇 Complexity and TPU’s Parallelism
Today’s parallel computing devices, including GPUs and TPUs, often consist a heterogeneous set of pipelined compute units U = {𝑈 1, 𝑈 2, . . . , 𝑈𝐾 }, where each 𝑈𝑘 denotes a class of units. Let 𝑃𝑘 be the number of parallel instances of unit type 𝑈𝑘 . We decompose the total work of size 𝑁 into per-unit contributions 𝑊𝑘 (𝑁 ), 𝑘 ∈ [1, . . . , 𝐾], where 𝑊𝑘 (𝑁 ) counts the operations mapped to unit 𝑈𝑘 . All these on-chip compute units are deeply pipelined, the execution time is bounded by the slowest (bottleneck) pipeline stage. We define the BIG-𝑇 complexity (termed as Big-𝑇 notation) of an algorithm on a parallel computing device as 𝑊 𝑇 (𝑁 ) = 𝑂 max max 𝑃 𝑘 , 𝑀𝑒𝑚 𝑘
𝑘
if there exists anconstant 𝑐 > 0 and o 𝑁 0 such that for all 𝑁 ≥ 𝑁 0 , 𝑊 𝑊 𝑇 (𝑁 ) ≤ 𝑐 · max max𝑘 𝑃 𝑘 , 𝑀𝑒𝑚 . Here, the term “max𝑘 𝑃 𝑘 ” cap𝑘 𝑘 tures the bottleneck among heterogeneous pipelined compute units. “𝑀𝑒𝑚” captures the overall latency of off-chip data access. Together, they provide an asymptotic lower latency bound of a proposed algorithm on parallel computing devices. Under sufficiently large workload, the ideal parallelism of compute units in a TPU used for Big-𝑇 analysis is listed in Fig. 2.
3.2
Arithmetic Optimizations
Table 1: Big-𝑇 Complexity of Arithmetic (Bottleneck in Red) Kernel Radix Mont. MXU RNS Lazy
VPU Span 𝐷2 𝑃𝐴𝑅VPU 4𝐷 𝑃𝐴𝑅VPU
MXU Span 𝐷2 𝑃𝐴𝑅MXU 𝐷2 𝑃𝐴𝑅MXU
XLU Span 𝐷 2 log 𝐷 𝑃𝐴𝑅S ∼0
Memory Span 𝐷 𝐵𝑊HBM 2𝐷 𝐵𝑊HBM
3.2.1 Challenge: Radix Montgomery Reduction is Dominated by Memory Reorganization. For a value with 𝐷 32-bit digits, radix-232 Montgomery multiplication requires (1) a 𝑂 (𝐷 2 ) big-integer multiplication and (2) a 𝑂 (𝐷 2 ) Montgomery reduction, each containing a sequential 𝐷-step carry-propagation chain. Carry propagation demands fine-grained shuffling and index permutation on a TPU, 3
Actual Value Range Changes In Actual range of Value Finite Field in ZKP
Computation Illustration 32D-bit x
MXU-centric ModMul
Execution Order
Extended Finite Field
Mul. Result
Partially Reduced
x y
xy MXU-centric Lazy Reduction our contribution xy
MXU-centric RNS ModMul Overall Result result Final Full Modular Reduction Final Reduction result <
Algorithm 1 MXU-Centric RNS Lazy Reduction (MXU RNS Lazy)
32D-bit y
[𝑥]𝑏 , 𝑏 ∈ [0, 𝐵). ByteDecompose(𝑥)→ 1: 2: for 𝑏 = 0 to 𝐵 − 1: 3: 𝑥𝑏 ← 𝑥 [8𝑏 : 8(𝑏 + 1)] ⊲ Select b-th byte of 𝑥, full parallel 4: Return [𝑥]𝑏 , 0 ≤ 𝑏 < 𝐵 ByteMerge([𝑥]ℎ , ℎ ∈ [0, 𝐻 )) → 𝑥 5: for ℎ = 0 to 𝐻 − 1: 6: 𝑥 [8ℎ : 8(ℎ + 1)] = [𝑥]ℎ ⊲ Parallel among H bytes 7: Return 𝑥 Precomputation(𝑄, 𝑃, 𝑤) Require: 𝑄; 𝑃; 𝑤: the Montgomery factor used in element wise reduction; 𝑦 = (2𝑤 ) −1 , 𝑧 = 2𝑤 , 𝑦 and 𝑧 are co-prime. 8: Initialize 𝐼𝑖 , 𝐼𝑖,𝑏 , [𝐸] 𝑖,𝑏,𝑗,ℎ , [𝑓 ] 𝑖 . 9: 𝐼𝑖 = ( 𝑞𝑄𝑖 ) −1 mod 𝑝𝑖 ( 𝑞𝑄𝑖 𝑦) mod 𝑄 𝑄 𝑄 10: 𝐼𝑖,𝑏 = ( 𝑞 ) −1 mod 𝑝𝑖 ( 𝑞 𝑦) << 8𝑏 mod 𝑄 ⊲ BAT[35] 𝑖 𝑖 11: [𝐸] 𝑖,𝑏,𝑗,ℎ = ByteDecompose 𝑧(𝐼𝑖,𝑏 mod 𝑀) mod 𝑝 𝑗 ℎ l 𝑢m 𝐼 2 12: [𝑓 ] 𝑖 = 𝑖𝑄 13: [𝑔] 𝑗 = (−𝑧(𝑄 mod 𝑀)) mod 𝑝 𝑗 14: Return [𝐸] 𝑖,𝑘 , [𝑓 ] 𝑖 , [𝑔] 𝑗
in
32-bit VecModMul in
8-bit MatMul +
32-bit VecOp in
RNS representation of value value
in
2D 32-bit moduli
Figure 3: Illustration of actual value change and computational flow of MXU-centric RNS Lazy Modular Multiplication. which appears as the XLU span bottleneck in Tab. 1. Thus, the limiting factor is not compute throughput, but the repeated memory reorganization inherent in Radix Montgomery Reduction. 3.2.2 Solution: Reducing 𝑂 (𝐷 2 ) → 𝑂 (2𝐷) for Big-Integer Multiplication via RNS. To eliminate quadratic-cost digit multiplications, we encode each high-precision prime-field element in an extended non-prime field F𝑄 following [19]. We select 𝑄 > 𝑀 2 to guarantee that for any 𝑥, 𝑦 ∈ F𝑄 , the product of their RNS vectors never overflows 𝑄, avoiding a true reduction by 𝑄. Hence, the multiplication is fully decomposed into independent 32-bit limb multiplications. This transformation converts the compute pattern from quadratic all-to-all digit multiplications into linear digit-independent multiplications.
Main(𝑥𝑄 , 𝑄, 𝑃, w) ⊲ MXU-centric RNS Lazy Reduction Î Î Require: 𝑄 = 𝑖 𝑞𝑖 ; 𝑃 = 𝑗 𝑝 𝑗 ,𝑖 ∈ [0, 𝐼 ), 𝑗 ∈ [0, 𝐽 ); w; pick Í 𝑢 ≥ ⌈𝑙𝑜𝑔2( 𝑖 𝑞𝑖 )⌉. 𝑥𝑄 = CRT𝑄 (𝑥) is RNS representation of 𝑥 in F𝑄 , containing 𝐼 residues with 𝐵 Bytes each. When being written as 𝑥𝑄 𝑖 , each iteration specifies entire 32-bit residue value of 𝑥𝑄 . When being written as 𝑥𝑄 𝑖,𝑏 , 𝑏 ∈ [0, 𝐵), each iteration specify 𝑏-th byte of 𝑖-th residue. [𝑥 𝑃 ] 𝑗,ℎ has 𝐽 residues with 𝐻 bytes each. Note: we refer to [24] for details of CRT𝑄 . Ensure: 𝑥 𝑃 = ((((𝑥 × 𝑦) mod 𝑀) mod 𝑄) × 𝑧) mod 𝑃 in F𝑃 15: [𝐸] 𝑖,𝑏,𝑗,ℎ , [𝑓 ] 𝑖 , [𝑔] 𝑗 ← Precomputation(𝑄, 𝑃, 𝑤) 16: 𝑣 ← 𝑥𝑄 · [𝑓 ] 𝑖 ⊲ Dot Product, fused into L18 as one MatMul 𝑣 𝑖 17: 𝑘 ← 2𝑢 ⊲ 32-bit Vectorized Shift 18: 𝑥𝑄 = ByteDecompose( 𝑥 ) 𝑄 𝑖 𝑖,𝑏 Í 19: [𝑥 𝑃 ] 𝑗,ℎ = 𝑥𝑄 𝑖,𝑏 × [𝐸] 𝑖,𝑏,𝑗,ℎ ⊲ (Einsum) uint8 MatMul 20: [𝑥 𝑃 ] 𝑗 = ByteMerge([𝑥 𝑃 ] 𝑗,ℎ ) 21: Return [𝑥 𝑃 ] 𝑗 + 𝑘 · [𝑔] 𝑗 ⊲ 32-bit ScalarVecMul + VecAdd
3.2.3 Solution: Reducing 𝑂 (𝐷 2 ) Reduction via MXU-Centric RNS Lazy Reduction. Although multiplication is linearized in F𝑄 , modular reduction by 𝑀 must still be performed over the prime field F𝑀 . A naïve approach would require RNS reconstruction from F𝑄 → prime-field reduction (mod 𝑀) → RNS decomposition into F𝑄 . We reduce this overhead by applying Basis Aligned Transformation [35] to Simon’s reduction [24], collapsing this multi-stage pipeline into one 8-bit matrix multiplication (Lines 15 and 18 in Alg.1) plus four 32-bit vector operations (Lines 16,17,19,20 in Alg.1). This shifts the dominant 𝑂 (𝐷 2 ) term to a low-precision matrix multiplication, with its cost amortized by MXU, yielding 𝑂 (𝐷) latency overhead. It also removes all carry-propagation from the critical path, making it throughput-bound rather than shuffle-bound. 3.2.4 Walk-Through Example. Given 𝑥, 𝑦 ∈ F𝑀 , we convert them into extended-field RNS form 𝑥𝑄 , 𝑦𝑄 , each containing 2𝐷 32-bit digits. The modular multiplication proceeds as: • RNS-level 32-bit VecMul with Montgomery reduction, 𝑂 (2𝐷). • MXU-centric RNS Lazy Reduction (Lines 14–20 of Alg.1, 𝑂 (𝐷 2 )). All sequential carry propagation are removed. Computation becomes dominated by vectorized operations, as highlighted in Tab. 1. In summary, MORPH encapsulates computation over F𝑀 inside a larger field F𝑄 , paying (1) a 2× memory footprint for the extended RNS representation, and (2) additional RNS reconstruction/decomposition. These costs are amortized among sea-of-MACs in MXU, such that the bottleneck shifts from sequential shufflebound carry propagation to VPU-bounded vectorized operations, enabling better throughput/latency than Radix Mont.
3.3
and (3) memory overhead from loading duplicated points of different windows from off-chip memory. MORPH introduces Layout-Stationary Pippenger (LS-PPG), which enforces a common sharding and layout across presorting and Bucket Accumulation (BA), the latency-dominant phase of MSM. LSPPG shards reduction-free dimensions across tensor cores, allowing each core to process an independent partition. This makes presorting a direct producer of BA-ready points, eliminating cross-TPU communication, intra-TPU layout reorganization. By fusing presorting with BA, post-sorted points are consumed immediately upon generation rather than being written back to and later reloaded from off-chip HBM. Beyond BA, MORPH performs Bucket Reduction in a tree form to expose more parallelism and improve VPU compute utilization. Alg. 2 and Fig. 4 illustrate LS-PPG. Consequently, the ideal memory span is reduced from 𝐵𝑊𝐾𝐻𝑁𝐵𝑀 to a single pass over scalars and points, i.e., 𝐵𝑊2𝑁 in Tab. 2. 𝐻 𝐵𝑀
Dataflow Optimizations
3.3.1 Multi-Scalar Multiplication (MSM). Challenge. porting SotA MSM to TPU leads to three inefficiencies: (1) inter-TPU communication, (2) intra-TPU layout reorganization, 4
S
-th Scalar;
P
-th Point;
2 bits out of 256 bits The 2nd Window c=2 (One window) S0 P0 P1 S1 SN-1
-th Bucket in -th window;
11 P0 01 P1 10 PN-1
Windows
Parallelism: Layout Transformation:
0 (00)
P3,P8,P11
1 (01) 2 (10)
P1,P4,P5 P9,P7,PN-1
3 (11)
P0,P7,P10
W
Accumulated Result of -th Window
Bucket Accumulation (BA)
Bucket Reduction (BR)
Batch Point Add
Accumulated Result Per Window
Buckets(Value) Points in bucket
PN-1
256-bit
B , Pre sorting
0 P3 1 P1
P8 P4
P11 P5
B0,2
2 P9 3 P0
P7 P7
PN-1 P10
B2,2 B3,2
K Points
2c
2c
MSM Result
drop zero bucket
B1,2
S0 S1
P0 P1
SN-1
K(2c-1)/c None
( )
Scatter
Window Merge (WM)
2c
None
PN-1
1(
)
Reshape (TPU)
Figure 4: Illustration of LS-PPG’s (Alg. 2). Takeaway: MORPH minimizes data re-sharding and layout reorganization.
N=RC
ω2 ω4 ω6 ω8
ω4 ω8
ω8 ω8 ω8
P:
R R
TF
C
×R
R
Matrix ModMul (R,R,C)
N/2 N/2 N/2 Parallelism: R*N LT: Shuffle Shuffle Shuffle Layout: Reshape -> 2D
(a) Butterfly NTT
C
TF
C
× C TF
R
Vector ModMul (R,C)
R=R1 R2 R1
C R1
Matrix ModMul (R,C,C)
N
N*C
None
Flatten -> 1D
TF
R2
R2
R2 R1 TF
× R1
C
C
ω8
C
ω4 ω8
R1
Matrix ModMul Vector ModMul batch C - (R1,R1,R2) batch C - (R1,R2)
R2
× R2 TF
R
C R
Matrix ModMul batch C - (R1,R2,R2)
TF
Vector ModMul (R,C)
C R
C
× C TF Matrix ModMul (R,C,C)
Parallelism: R1*N
N
R2*N
N
N*C
Layout: Reshape -> 3D
None
None
Reshape -> 2D
Flatten -> 1D
(b) 3-step NTT
(c) 5-step NTT (Proposed by This Work)
Figure 5: Performance Analysis of Different NTT. (𝑁 = 𝑅 × 𝐶, 𝑅 = 𝑅1 × 𝑅2 ), TF refers to twiddle factors, which are generated offline. Takeaway: MORPH further recursively reduces the row-wise NTT in (b) into an 3-step NTT to reduce overall computation. Table 2: Big-𝑇 Complexity of Dataflows (Bottleneck in Red)
Algorithm 2 Layout Stationary Pippenger (LS-PPG)
VPU Span and MXU Span
XLU Span
Memory Span
LS-PPG
𝐾𝑁 2𝐾 (2𝑐 − 1) (𝐾 − 1)(1 + 𝐶) + 1 + + 𝑃𝐴𝑅 𝐵𝐴 𝑃𝐴𝑅 𝐵𝑅 𝑃𝐴𝑅𝑊 𝑀 𝐾𝑁 4𝐾 (2𝑐 − 1) (𝐾 − 1) (1 + 𝐶) + 1 + + 𝑃𝐴𝑅 𝐵𝐴 𝑃𝐴𝑅 𝐵𝑅_𝑛𝑒𝑤 𝑃𝐴𝑅𝑊 𝑀
2𝑐 · 𝑁 log 𝑁 𝑃𝐴𝑅S 2𝑐 · 𝑁 log 𝑁 𝑃𝐴𝑅S
𝐾 ·𝑁 𝐵𝑊HBM 2𝑁 𝐵𝑊HBM
NTT
VPU Span
XLU Span
Memory Span
MSM
Require: Initialize 𝑀𝑆𝑀, 𝑆𝑛 , 𝑃𝑛 ,𝑊𝑘 , 𝐵𝑘,𝑗 ,𝑊𝑘 as 0, 𝑗 ∈ [0, 2𝑐 − 1], 𝑘 ∈ [0, 𝐾], 𝑁 ′ (max point count per 𝐵𝑘,𝑗 ), definition in Fig. 4. Bucketize 1: sharding for 𝑘 ∈ 𝐾 ⊲Distributed across multi-TPUs 2: Initialize 𝑃𝑘′ ← O of shape [2𝑐 , 𝑁 ′ ] ⊲ 2𝑐 buckets. 3: Π𝑘 ← argsort(𝐼𝑘 ) ⊲𝐼𝑘 : 𝑘-th slice of all scalars 𝑆𝑛 Gather 𝑃𝑛 based on Π𝑘 , and scatter them into 𝐼𝑘 [𝑛]-th bucket 4: in 𝑘-th window (𝑊𝑘 ). Bucket Accumulation (BA) 5: Initialize bucket 𝐵𝑘,𝑗 with zero points. 6: sharding for 𝑘 ∈ 𝐾 ⊲Distribute windows across multi-TPUs 7: parallel for 𝑗 ∈ [1, 2𝑐 − 1] 8: for 𝑛 ′ ∈ 𝑁 ′ ′ 9: 𝐵𝑘,𝑗 += 𝑃𝑘,𝑛 ′ ,𝑗 Bucket Reduction (BR) ⊲ Tree-based Accumulation 10: Initialize 𝑊𝑘,𝑗 ← O for all 𝑘 ∈ [0, 𝐾 − 1] and 𝑗 ∈ [0, 2𝑐 − 1] 11: for 𝑠 = 𝑐 − 1 down to 0 ⊲level of tree (𝐿) (𝑅) 12: Let 𝐵𝑘,𝑏 ← 𝐵𝑘,2𝑏 and 𝐵𝑘,𝑏 ← 𝐵𝑘,2𝑏+1 for 𝑏 ∈ [0, 2𝑠 − 1]
Presort-PPG
MXU Span
𝑁 log 𝑁 𝑁 log 𝑁 (𝑁 + 𝑁 ) 0 𝑃𝐴𝑅VPU 𝑃𝐴𝑅S 𝐵𝑊HBM 𝑁 𝑁 (𝑅 + 𝐶) 2𝑁 (2𝑁 + 𝑅 2 + 𝐶 2 ) 3-step NTT 𝑃𝐴𝑅VPU 𝑃𝐴𝑅MXU 𝑃𝐴𝑅T 𝐵𝑊HBM (2𝑁 + 𝑅12 + 𝑅22 + 𝑅 + 𝐶 2 ) 2𝑁 𝑁 (𝑅1 + 𝑅2 + 𝐶) 3𝑁 5-step NTT 𝑃𝐴𝑅VPU 𝑃𝐴𝑅MXU 𝑃𝐴𝑅T 𝐵𝑊HBM Note: 𝑃𝐴𝑅 𝐵𝐴 /𝑃𝐴𝑅 𝐵𝑅 /𝑃𝐴𝑅𝑊 𝑀 refers to 𝑃𝐴𝑅𝑉 𝑃𝑈 or 𝑃𝐴𝑅𝑀𝑋𝑈 in BA,BR,WM phase. 𝑃𝐴𝑅 𝐵𝑅_𝑛𝑒𝑤 =𝑐 . 𝑃𝐴𝑅 𝐵𝑅 =2. VPU / MXU runs the same number of PADD. Butterfly
3.3.2 Number-Theoretic Transformation (NTT). Challenge of Butterfly NTT and 3-step NTT. Tab. 2 shows the Big-𝑇 complexity of NTT variants. Running Butterfly NTT on a 3 VPU TPU shows 𝑃𝐴𝑅 𝑃𝐴𝑅S ≈ 𝑂 (10 ), making XLU the critical bottleneck slowing down butterfly NTT by 𝑂 (103 ). This is slow even though it offers minimal compute complexity. The 3-step NTT [35] removes shuffles by reforming a 𝑁 -degree NTT into two (𝑅, 𝑅, 𝐶) matrix multiplications and 𝑁 -length vector operations, where 𝑁 = 𝑅 · 𝐶. However, its Big-𝑇 is dominated by the MXU span 𝑁 (𝑅 + 𝐶)/𝑃𝐴𝑅MXU , which grows with √ 𝑁 by a factor of (𝑅 + 𝐶). Even with balanced factors (𝑅 = 𝐶 = 𝑁 ), it incurs a memory span of at least 4𝑁 /𝐵𝑊HBM . Solution: 5-step NTT to reduce MXU span while preserving MXU utilization. To further reduce the MXU and memory span, MORPH introduces a 5-step NTT that replaces the row-wise 𝑅-degree NTT (step 1 in the 3-step NTT of Fig. 5b) with a 3-step NTT over (𝑅1, 𝑅2, 𝐶), as shown in Fig. 5c, where 𝑅 = 𝑅1 · 𝑅2 . This decomposes the large GEMM into multiple smaller GEMMs, reducing the effective MXU span by a factor of 𝑁 1/4 while being sufficiently large to fully utilize the MXU. The 5-step NTT is detailed in Eq. 1. h h i i 𝐶 ·𝑅1 𝑅 𝐶 2 𝑇 𝐹𝑅×𝑅 @ 𝑇 𝐹𝑅𝐶2×𝑅 ×𝑅2 @ 𝑎𝑅1 ×𝐶 ×𝑅2 @𝑇 𝐹𝑅1 ×𝑅1 ·𝑇 𝐹𝑅1 ×𝑅2 ·𝑇 𝐹𝐶 ×𝐶 (1)
(𝐿) (𝑅) Let 𝑊𝑘,𝑏 ←𝑊𝑘,2𝑏 and 𝑊𝑘,𝑏 ←𝑊𝑘,2𝑏+1 for 𝑏 ∈ [0, 2𝑠 − 1] 𝑠 14: sharding for 𝑏 = 0 to 2 − 1 15: parallel for 𝑘 = 0 to 𝐾 − 1 (𝐿) (𝑅) (𝑅) 16: 𝑊𝑘,𝑏 ← 𝑊𝑘,𝑏 + 𝑊𝑘,𝑏 + 𝐵𝑘,𝑏 (𝐿) (𝑅) 17: 𝐵𝑘,𝑏 ← 2 · 𝐵𝑘,𝑏 + 𝐵𝑘,𝑏 18: 𝑊𝑘 ← 𝑊𝑘,0 ⊲ root of the tree Window Merging (WM) 19: for 𝑘 in (𝐾, 0, -1): ⊲Reverse Window Index 𝑘 20: 𝑀𝑆𝑀 ×= 2𝑐 21: 𝑀𝑆𝑀 += 𝑊𝑘 22: Return 𝑀𝑆𝑀. Note 1: Coordinate and bytes dimensions are hidden for simplicity. Note 2: O is the zero point (0, 1, 1, 0) in Twisted Edward curves. (𝐿) (𝑅) (𝐿) (𝑅) ,𝑊𝑘,𝑏 ,𝑊𝑘,𝑏 are introduced for tree reduction. Note 3: 𝐵𝑘,𝑏 , 𝐵𝑘,𝑏 Note 4: N’ is the maximal number of points among all buckets.
13:
Here, @ denotes matrix multiplication and · represents element𝑘 wise multiplication. 𝑇 𝐹𝑅×𝐶 is a matrix [((𝜔𝑛 )𝑘 )𝑖 𝑗 ], 𝑖 ∈ [0, 𝑅), 𝑗 ∈ [0, 𝐶). 𝜔𝑛 : primitive 𝑁 -th root of unity. 𝑎 is the input tensor. 5
102 101
NTT
107 105 103
Speedup
Table 3: Latency Comparison (MORPH -tpuv6e8 v.s. GZKPV100). Unit: 𝜇s/ms for NTT/MSM. Improvement in red.
102
Workload
101
100 100 218 219 220 100 214 216 218 220 222 MSM Degree (TPUv6e8) NTT Degree (TPUv6e) Figure 6: MORPH Ablation Study under different degrees.
NTT (𝜇s)
MSM (ms)
Workload: We take the most popular primitives from zk-SNARKs (degree usually ranges from 214 ∼ 226 ) [3, 9, 15, 34]. We adopt the same evaluation setup as GZKP for MSM and NTT [1, 27], where we assume data comes in affine representation [2] and offline converted into Twisted Edwards form to minimize compute overhead. TPU Baseline: We adopt SotA GPU algorithms as our baseline, including (1) radix Montgomery Reduction [20, 27], (2) Presorting Pippenger [20], and (3) 3-step NTT [35], detailed in §2.3.3. MORPH TPU Setup: MORPH uses (1) MXU-centric RNS Lazy Reduction, (2) LS-PPG, and (3) 3-step and 5-step NTT. MORPH is implemented in JAX [10] and captures wall-clock latencies using the JAX Profiler (XProf) [5]. We use Google TPUv6e as platforms, and scale the number of chips to match V100’s power consumption.
𝑁 = 222 𝑁 = 224 𝑁 = 226
Throughput ↑
Latency (ms, log scale)
100 10 1 0 0
753-bit GZKP MORPH 150 15.1 490 52.8 910 337 7,460 2,195 33,670 13,107 141,400 OOM 2.6–10× 2.66 2.24 11.30 8.91 40.70 35.64 1.14–1.2×
Precision
256-bit 377-bit 753-bit
Algorithm
Radix Mont. MXU RNS Lazy 20 22 24 26 28 210 212 214 216
256-bit GZKP MORPH 50 11.6 90 24.8 280 109 1,070 716 4,960 4,694 20,990 21,382 0.98–4.3× 0.24 1.61 1.10 6.42 4.00 25.64 0.15–0.16
1k
V5p (753 bits) V6e (753 bits) V5p (256 bits) V6e (256 bits) 100 2 4 8 16 32 64 128 256 512
Batch Size (NTT) Batch Size (ModMul) Figure 7: ModMul and NTT under different batch sizes.
carry chains, whereas TPU latency grows only 1.3× ∼ 3× because RNS arithmetic maps efficiently onto the MXU. TPUs also avoid butterfly shuffling entirely. Their cost is determined by the compute spans in Tab. 2, including MatMul (𝑁 (𝑅+𝐶))/𝑃𝐴𝑅MXU , XLU twiddle steps (2𝑁 ∼ 3𝑁 ), and (2𝑁 + 𝑅 2 + 𝐶 2 )/𝐵𝑊HBM memory span, all of which are small for lower-degree NTTs.
MORPH Evaluation
We conduct three-stage ablation for MSM and NTT (Baseline, +Arithmetic, +Ari. + Dataflow, ) with latency and speedup plotted in Fig. 6. 4.2.1 Arithmetic Optimization Evaluation. Arithmetic optimizations yield 50 ∼ 90× latency reduction for MSM, collapsing BR/WM to near-constant cost. These gains come from MXU-centric RNS Lazy Reduction, which removes sequential carry propagation and reformulating modular reduction of big integer arithmetic into lowprecision dense matrix multiplication, allowing the MXU to amortize native 𝑂 (𝐷 2 ) computations into 𝑂 (𝐷) latency. Although this increases memory footprint by up to 2×, MSM remains computebound on VPU operations, so off-chip latency is effectively hidden. Overall, MXU-accelerated RNS-lazy arithmetic is the primary contributor to performance improvement for both MSM and NTT. 4.2.2 MSM – Dataflow Optimization Evaluation. MORPH enables up-to 3.1× speedup. Across all degrees, BA (Bucketized included) dominates the runtime, because of random scatter and gather. Dataflow optimizations reduce BR by ∼ 6× from parallelism increase. 4.2.3 NTT – Dataflow Optimization Evaluation. 3-step NTT is faster at lower degrees, while 5-step becomes faster at higher degrees with up-to 1.07× speedup. 5-step NTT forces the overall degree 𝑁 to be decomposed into three factors. At small degree, these factors are typically ≤ 26 = 64, leading to poor utilization of the (8, 128) VReg shape and reducing utilization for the MXU and VPU. As the degree grows (> 222 ), these factors become large enough to fully utilize VRegs, and 5-step benefits from reduced storage overhead of twiddle parameters.
4.3
214 216 218 220 222 224
Throughput ↑
4 Experiments 4.1 Experiment Setup
4.2
Scale
Latency (ms)
103
Speedup (x)
Speedup
Baseline +Arithmetic Opt +Ari.+DataFlow Opt
WM
Speedup (x)
102
BR
Runtime (ms, log scale)
104
BA Baseline +Arithmetic +Ari.+DataFlow Opt
Runtime (ms, log scale)
106
4.3.2 MSM. Tab. 3 reports estimated latency relative to GZKP. MORPH ’s LS-PPG allows TPUv6e8 to achieve competitive throughput for 753-bit MSM. The latency is dominated by Bucketize, whose per-point scatter exhibits poor memory-bandwidth utilization. This overhead becomes more severe with finer-grained point scatter, leading to lower performance than GZKP at 256-bit. A dedicated hardware reordering unit that hides fine-grained scatters behind computation could mitigate this bandwidth bottleneck [36, 37].
4.4
Ablation Study - Batch Size
MXU-centric RNS-Lazy sustains 4∼157× lower latency than Radix Montgomery across all batch sizes and precisions, and the performance gap widens as precision increases (256→377→753 bits) and as batch size increases, it further amplifies MXU utilization (Fig. 7). Increasing batch size from 1 to 128 gives 3.8/3.2× and 5.4/3.9× speedup on NTT (753 and 256 bit) for TPUv5p/v6e (Fig. 7). NTT favors≥ 128 batches to fully utilize VRegs. Latency plateaus beyond 128 batch size as MXU/VPU achieves its peak compute utilization.
5
Conclusion
This paper presents MORPH, the first framework to accelerate ZKP kernels on AI ASICs, deployed on Google TPU. Using a BigT formulation, we identify the two core bottlenecks and show how to overcome them: bridging high-precision modular arithmetic (>256-bit) to low-precision matrix/vector engines via an extended non-prime RNS field, and eliminating or hiding layoutreorganization and inter-device communication overhead through dataflow with layout-invariant sharding over non-reduction dimensions. MORPH makes TPU the SotA commodity platform for high-precision NTT throughput and position AI ASICs as a practical foundation for full ZKP protocol acceleration.
MORPH vs. SotAs
4.3.1 NTT. Tab. 3 highlights two core advantages of TPUs for largeprecision NTTs. First, TPUv6e8 achieves SoTA throughput in NTT: Its leading energy efficiency transfers to 10× higher throughput than GZKP (V100) for 753-bit NTTs. Second, TPUs scale more gracefully with precision. Increasing the modulus from 256- to 753-bits raises GPU latency by 6× ∼ 7× due to long Montgomery 6
5.1
Acknowledgments
[23] Kim Laine, Rachel Player, and Hao Chen. 2018. Microsoft SEAL: A Homomorphic Encryption Library. In Proceedings of the IEEE Symposium on Security and Privacy Workshops (SPW). IEEE, 123–126. [24] Simon Langowski and Srinivas Devadas. 2025. Efficient Modular Multiplication Using Vector Instructions on Commodity Hardware. Cryptology ePrint Archive, Paper 2025/1068. https://eprint.iacr.org/2025/1068 [25] Changxu Liu et al. 2024. Gypsophila: A Scalable and Bandwidth-Optimized Multi-Scalar Multiplication Architecture. In Proceedings of the 61st ACM/IEEE Design Automation Conference (San Francisco, CA, USA) (DAC ’24). Association for Computing Machinery, New York, NY, USA. doi:10.1145/3649329.3658259 [26] Changxu Liu, Hao Zhou, Patrick Dai, Li Shang, and Fan Yang. 2023. PriorMSM: An Efficient Acceleration Architecture for Multi-Scalar Multiplication. In Proceedings of the ACM/SIGDA International Symposium on Field-Programmable Gate Arrays. doi:10.1145/3678006 https://dl.acm.org/doi/10.1145/3678006. [27] Tianyu Ma, Zhen Zhang, Yuhao Zhang, and G. Edward Suh. 2023. gZKP: GPUAccelerated Zero-Knowledge Proof Generation. In Proceedings of the IEEE International Symposium on High-Performance Computer Architecture (HPCA). IEEE. [28] Weiliang Ma, Qian Xiong, Xuanhua Shi, Xiaosong Ma, Hai Jin, Haozhao Kuang, Mingyu Gao, Ye Zhang, Haichen Shen, and Weifang Hu. 2023. GZKP: A GPU Accelerated Zero-Knowledge Proof System. In Proceedings of the 28th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 2 (ASPLOS 2023). Association for Computing Machinery, Vancouver, BC, Canada, 340–353. doi:10.1145/3575693.3575711 [29] Ayumi Ohno, Kotaro Shimamura, and Shinya Takamaeda-Yamazaki. 2025. Accelerating Elliptic Curve Point Additions on Versal AI Engine for Multi-scalar Multiplication. arXiv:2502.11660 [cs.AR] https://arxiv.org/abs/2502.11660 [30] Xander Pottier, Thomas de Ruijter, Jonas Bertels, Wouter Legiest, Michiel Van Beirendonck, and Ingrid Verbauwhede. 2025. OPTIMSM: FPGA hardware accelerator for Zero-Knowledge MSM. IACR Transactions on Cryptographic Hardware and Embedded Systems 2025, 2 (2025), 489–510. [31] Andy Ray, Benjamin Devlin, Fu Yong Quah, and Rahul Yesantharao. 2023. Hardcaml MSM: A High-Performance Split CPU-FPGA Multi-Scalar Multiplication Engine. In Proceedings of the ACM Symposium on Field-Programmable Gate Arrays. doi:10.1145/3626202.3637577 https://dl.acm.org/doi/10.1145/3626202.3637577. [32] Brandon Reagen, Woojoo Choi, David Brooks, Gu-Yeon Wei, and Hsien-Hsin S. Lee. 2021. HEAX: An Architecture for Computing on Encrypted Data. In Proceedings of the ACM/IEEE International Symposium on Computer Architecture (ISCA). IEEE, 1113–1126. [33] Nikola Samardzic, Simon Langowski, Srinivas Devadas, and Daniel Sanchez. 2024. Accelerating Zero-Knowledge Proofs Through Hardware-Algorithm Co-Design. In 2024 57th IEEE/ACM International Symposium on Microarchitecture (MICRO). 366–379. doi:10.1109/MICRO61859.2024.00035 [34] Roman Storm, Alexey Pertsev, and Roman Semenov. 2020. Tornado Cash Privacy Solution. https://tornado.cash/Tornado.cash_whitepaper_v1.4.pdf White Paper. [35] Jianming Tong, Tianhao Huang, Jingtian Dang, Leo de Castro, Anirudh Itagi, Anupam Golder, Asra Ali, Jeremy Kun, Jevin Jiang, Arvind, G. Edward Suh, and Tushar Krishna. 2026. Leveraging ASIC AI Chips for Homomorphic Encryption. In 2026 IEEE International Symposium on High Performance Computer Architecture (HPCA). 1–18. doi:10.1109/HPCA68181.2026.11408507 [36] Jianming Tong, Anirudh Itagi, Parsanth Chatarasi, and Tushar Krishna. 2024. FEATHER: A Reconfigurable Accelerator with Data Reordering Support for LowCost On-Chip Dataflow Switching. In Proceedings of the 51th Annual International Symposium on Computer Architecture (Argentina) (ISCA ’24). Association for Computing Machinery, Argentina. [37] Jianming Tong, Yujie Li, Devansh Jain, Charith Mendis, and Tushar Krishna. 2026. MINISA: Minimal Instruction Set Architecture for Next-gen Reconfigurable Inference Accelerator. In Proceedings of the 34th Annual International Symposium on Performance Analysis of Systems and Software (Seoul, Korea) (ISPASS ’26). [38] Cheng Wang and Mingyu Gao. 2025. UniZK: Accelerating Zero-Knowledge Proof with Unified Hardware and Flexible Kernel Mapping. In Proceedings of the 30th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 1 (Rotterdam, Netherlands) (ASPLOS ’25). Association for Computing Machinery, New York, NY, USA, 1101–1117. doi:10.1145/3669940.3707228 [39] Samuel Williams, Andrew Waterman, and David Patterson. 2009. Roofline: an insightful visual performance model for multicore architectures. Commun. ACM 52, 4 (April 2009), 65–76. doi:10.1145/1498765.1498785 [40] Zhengbang Yang et al. 2025. LegoZK: A Dynamically Reconfigurable Accelerator for Zero-Knowledge Proof. In 2025 IEEE International Symposium on High Performance Computer Architecture (HPCA). doi:10.1109/HPCA61900.2025.00020 [41] Yancheng Zhang et al. 2025. zkVC: Fast Zero-Knowledge Proof for Private and Verifiable Computing. In Proceedings of the 62nd Annual ACM/IEEE Design Automation Conference (San Francisco, California, United States) (DAC ’25). IEEE Press. doi:10.1109/DAC63849.2025.11132681 [42] Ye Zhang, Shuo Wang, Xian Zhang, et al. 2021. PipeZK: Accelerating ZeroKnowledge Proof with a Pipelined Architecture. In Proceedings of the 2021 ACM/IEEE 48th Annual International Symposium on Computer Architecture (ISCA).
This work was supported in part by ACE, one of the seven centers in JUMP 2.0, a Semiconductor Research Corporation (SRC) program sponsored by DARPA. We thank reviewers for their feedbacks.
References [1] [n. d.]. Benchmark harness for FPGA MSM implementations in the ZPRIZE competition. https://github.com/z-prize/prize-gpu-fpga-msm/tree/main/harness. Accessed: April 1, 2025. [2] [n. d.]. XYZZ coordinates for short Weierstrass curves. https://www.hyperelliptic. org/EFD/g1p/auto-shortw-xyzz.html. Accessed: April 9, 2025. [3] [n. d.]. ZK Rollup Architecture. Online. https://zksync.io/faq/tech.html#zkrollup-architecture Accessed: April 2025. [4] 2025. Amazon Trainium. https://aws.amazon.com/ai/machine-learning/ trainium/. Accessed: [Insert Date Accessed, e.g., April 9, 2025]. [5] 2025. Perfetto Trace Viewer. https://perfetto.dev/. Accessed: [Insert Date Accessed, e.g., April 9, 2025]. [6] 2025. yrrid GPU MSM Library. https://github.com/yrrid/combined-msm-gpu. Accessed: [Insert Date Accessed, e.g., April 9, 2025]. [7] Kaveh Aasaraai, Don Beaver, Emanuele Cesena, Rahul Maganti, Nicolas Stalder, and Javier Varela. 2022. CycloneMSM: FPGA Acceleration of Multi-Scalar Multiplication. Technical Report. IACR. https://eprint.iacr.org/2022/1396.pdf. [8] Jacob Austin, Sholto Douglas, Roy Frostig, Anselm Levskaya, Charlie Chen, Sharad Vikram, Federico Lebron, Peter Choy, Vinay Ramasesh, Albert Webson, and Reiner Pope. 2025. How to Scale Your Model. Online. (2025). Retrieved from https://jax-ml.github.io/scaling-book/. [9] Eli Ben-Sasson, Alessandro Chiesa, Christina Garman, Matthew Green, Ian Miers, Eran Tromer, and Madars Virza. 2014. Zerocash: Decentralized Anonymous Payments from Bitcoin. Technical Report. Zerocash Project. http://zerocashproject.org/media/pdf/zerocash-extended-20140518.pdf [10] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. 2018. JAX: composable transformations of Python+NumPy programs. http://github.com/jax-ml/jax [11] Jeff Buss et al. 2021. Intel HEXL: High-Performance Homomorphic Encryption Primitives. In Proceedings of the Workshop on Encrypted Computing & Applied Homomorphic Cryptography (WAHC). [12] Wonseok Choi, Jongmin Kim, and Jung Ho Ahn. 2025. Cheddar: A Swift Fully Homomorphic Encryption Library Designed for GPU Architectures. In Proceedings of the 31st ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 1 (USA) (ASPLOS ’26). Association for Computing Machinery. doi:10.1145/3760250.3762223 [13] Alhad Daftardar, Brandon Reagen, and Siddharth Garg. 2024. SZKP: A Scalable Accelerator Architecture for Zero-Knowledge Proofs. In Proceedings of the 2024 International Conference on Parallel Architectures and Compilation Techniques (Long Beach, CA, USA) (PACT ’24). Association for Computing Machinery, New York, NY, USA, 271–283. doi:10.1145/3656019.3676898 [14] Wei Dai and Berk Sunar. 2015. cuHE: A Homomorphic Encryption Accelerator Library. IACR Cryptology ePrint Archive 2015 (2015), 1043. [15] Dark Forest Team. [n. d.]. Announcing Dark Forest. Blog post. https://blog.zkga. me/announcing-darkforest Accessed: April 2025. [16] Shengyu Fan, Zhiwei Wang, Weizhi Xu, Rui Hou, Dan Meng, and Mingzhe Zhang. 2022. TensorFHE: Achieving Practical Computation on Encrypted Data Using GPGPU. arXiv:2212.14191 [cs.AR] https://arxiv.org/abs/2212.14191 [17] Ariel Gabizon, Zachary J. Williamson, and Oana Ciobotaru. 2019. Plonk: Permutations over Lagrange-bases for Oecumenical Noninteractive arguments of Knowledge. IACR Cryptol. ePrint Arch. 2019 (2019), 953. https://eprint.iacr.org/2019/953. [18] Jens Groth. 2016. On the Size of Pairing-based Non-interactive Arguments. In Advances in Cryptology - EUROCRYPT 2016 (Lecture Notes in Computer Science, Vol. 9666). Springer, 305–326. doi:10.1007/978-3-662-49896-5_11 [19] David Jacquemin, Ahmet Can Mert, and Sujoy Sinha Roy. 2022. Exploring RNS for Isogeny-based Cryptography. Cryptology ePrint Archive, Paper 2022/1289. [20] Zhuoran Ji et al. 2024. Accelerating Multi-Scalar Multiplication for Efficient Zero Knowledge Proofs with Multi-GPU Systems. In Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 3 (La Jolla, CA, USA) (ASPLOS ’24). Association for Computing Machinery, New York, NY, USA, 57–70. doi:10.1145/3620666.3651364 [21] Norman P. Jouppi et al. 2023. TPU v4: An Optically Reconfigurable Supercomputer for Machine Learning with Hardware Support for Embeddings. arXiv:2304.01433 [cs.AR] https://arxiv.org/abs/2304.01433 [22] Jongmin Kim, Sangpyo Kim, Jaewan Choi, Jaiyoung Park, Donghwan Kim, and Jung Ho Ahn. 2023. SHARP: A Short-Word Hierarchical Accelerator for Robust and Practical Fully Homomorphic Encryption. In Proceedings of the 50th Annual International Symposium on Computer Architecture (Orlando, FL, USA) (ISCA ’23). Association for Computing Machinery. doi:10.1145/3579371.3589053 7