arXiv:2606.19866v1 [cs.CR] 18 Jun 2026
Low-Cost Multi-Precision Systolic Arrays for Accelerating FHE NTTs on AI ASICs George Alexakis
Dimitrios Schoinianakis
Giorgos Dimitrakopoulos
Electrical and Computer Engineering Democritus University of Thrace, Xanthi, Greece
Nokia Bell Labs Athens, Greece
Electrical and Computer Engineering Democritus University of Thrace, Xanthi, Greece
Abstract—Fully Homomorphic Encryption (FHE) ensures robust data privacy but suffers from prohibitive computational overhead. Accelerating FHE on AI hardware like Tensor Processing Units (TPUs) is promising, yet fundamentally limited by a precision mismatch: TPUs are optimized for 8-bit arithmetic, whereas FHE and its critical parts such as the Number Theoretic Transform (NTT), demand high precision. Current approaches bridge this gap using matrix decomposition to execute NTT computations on low-precision matrix engines. However, reconstructing the full-precision results requires shift-and-add accumulation that does not match the dataflow of matrix multiplication. This forces offloading full-precision reconstruction from matrix engines to vector processors that disrupts the matrix multiplication dataflow, creating significant performance bottleneck. To resolve this limitation, we propose a minimally modified multi-precision systolic array that performs full-precision output reconstruction natively within the array in sync with low-precision matrix multiplication under a uniform dataflow. Synthesized at 7nm with OpenRoad, our design incurs negligible hardware overhead. Cycle-accurate simulations using SCALE-Sim demonstrate that natively executing NTTs on the proposed architecture achieves at least 1.33× speedup, for transform sizes 212 to 216 on 128×128 matrix engines, successfully enabling standard AI hardware to support high-precision FHE acceleration.
I. I NTRODUCTION Fully Homomorphic Encryption (FHE) offers a powerful security guarantee: it enables arbitrary computation directly on encrypted data without ever exposing plaintext. Despite this promise, its extreme computational overhead has traditionally necessitated specialized, custom-built hardware [1], [2]. Designing such bespoke accelerators, however, is prohibitively expensive and difficult to scale. A compelling alternative is to leverage existing AI accelerators [3], [4], [5], such as Tensor Processing Units (TPUs), which already deliver massive parallelism and high energy efficiency [6]. Repurposing these platforms provides a practical path toward closing the performance gap that currently limits real-world FHE adoption [2]. Mapping FHE workloads to AI accelerators reveals a fundamental architectural mismatch: TPUs are optimized for 8-bit integer arithmetic [6], whereas FHE demands modular arithmetic over large integer moduli [7]. Although FHE mitigates this using the Residue Number System [8], which splits large integers into smaller, independent 30-60 bit channels, these components remain too large for native TPU matrix engines. To utilize TPU matrix engines, a common workaround is to decompose high-precision values into multiple low-precision (e.g., 8-bit) matrix operations, whose results must later be aggregated to full-precision. This reconstruction requires a
sequence of shift-and-accumulate operations and digit-level carry propagation. Current TPU architectures cannot perform these operations within their systolic matrix engines, forcing computation to spill over to external vector units. This “external accumulation” incurs substantial data movement overhead, disrupts matrix multiplication dataflow, and leaves the matrix engine underutilized, effectively negating the advantages of the accelerator. In this work, we show that efficient FHE execution on AI accelerators is achievable with minimal hardware changes. We propose a lightweight modification to the processing elements of the last row of a systolic array that enables in-place shiftand-add operations directly within the matrix engine. Unlike prior approaches for multi-precision matrix computation that rely on reconfigurable SIMD-style systolic arrays [9], our design introduces negligible hardware overhead while supporting multi-precision accumulation natively. By embedding reconstruction into the uniform dataflow of matrix multiplication, we eliminate the need for external accumulation and unlock sustained high utilization of the accelerator. This approach bridges the gap between the high-precision requirements of FHE and the low-precision efficiency of modern AI hardware, enabling practical deployment without sacrificing compatibility with existing workloads. The contributions of this work can be summarized as follows: We identify the full-precision reconstruction phase, specifically shift-and-add accumulation and digit-level carry propagation, as a critical performance bottleneck when mapping core FHE operations, such as the Number Theoretic Transform (NTT) used to accelerate polynomial multiplications, onto AI accelerators. • We introduce a minimally modified systolic array that performs full-precision reconstruction internally and in sync with low-precision matrix multiplication under a uniform dataflow. Supporting inside the systolic array this dual operation of low-precision matrix multiplication and full-precision result reconstruction adds negligible hardware cost (less than 1%), as quantified by implementation results at 7nm using OpenROAD [10]. • We evaluate performance using the SCALE-Sim cycleaccurate simulator. For the NTT, which dominates FHE runtime, we demonstrate at least 1.33× speedup for transform sizes ranging from 212 to 216 , while preserving compatibility with conventional low-precision AI matrix engines and incurring minimal overhead. •
II. BACKGROUND AND R ELATED W ORK FHE schemes are dominated by linear algebra and polynomial arithmetic [11]. Ciphertexts are represented as vectors of ring elements, and core operations such as key switching, automorphisms, and ciphertext multiplication ultimately reduce to structured linear transformations and polynomial convolutions [3], [4], [11], [5]. Recent work has shown that many of these kernels can be reformulated using dense matrix operations, enabling acceleration on GPUs and AI accelerators. TensorFHE [3] expresses FHE primitives for GPU tensor cores, while TensorFHE+ [4] extends this idea through a linear-algebra abstraction that improves portability across heterogeneous accelerators. CROSS [5] further demonstrates the feasibility of mapping FHE workloads onto systolic-array accelerators by reformulating high-precision modular arithmetic into lowprecision matrix multiplications. This transformation is particularly critical for core FHE operations [12], such as the NTT kernels that simplify polynomial multiplication. A. Matrix-Based NTT and Mapping to AI ASICs Custom FHE accelerators typically rely on dedicated parallel butterfly units to compute the NTT [2]. In contrast, approaches that target AI ASICs leverage existing matrix multiplication engines and vector processors. In this setting, the large NTTs required by FHE applications (ranging from 212 to 216 elements) are executed using the 4-step NTT algorithm [13], which decomposes the monolithic transform into a sequence of matrix-friendly sub-operations. To execute the 4-step NTT algorithm on an 1D input of size N = r × c the input is first reshaped into an r × c matrix. Then, the following four steps are executed: Column-wise Transforms: The algorithm performs c independent r-point NTTs along the columns of the matrix. This is executed as a matrix multiplication between a precomputed r × r transform matrix and the r × c data matrix. Twiddle Factor Multiplication: To correct for the phase interactions introduced by moving to two dimensions, the intermediate results undergo an element-wise multiplication with a precomputed r × c matrix of twiddle factors. Matrix Transpose: The r × c matrix is transposed into a c × r layout. This structural reorganization aligns the original row data into contiguous columns. Row-wise Transforms: Finally, r independent c-point NTTs are executed along the newly transposed columns. This is computed via a second matrix multiplication, this time using a c × c transform matrix. Ultimately, the 4-step replaces a single N √ decomposition √ point NTT with two N × N matrix multiplications, one element-wise multiplication, and one transpose operation. This formulation aligns well with the systolic matrix units found in modern AI accelerators: systolic arrays can natively execute the matrix multiplications, while auxiliary processing units handle the transpose and element-wise operations. B. Precision mismatch between FHE and AI ASICs Mapping FHE workloads onto AI ASICs exposes a significant architectural mismatch. FHE requires modular arithmetic
over large integer moduli (often spanning hundreds of bits), whereas modern AI ASICs rely on 8-bit integer units to maximize energy efficiency and scalability (typically 128×128 processing elements or larger) [6]. To bridge this precision gap, existing approaches, such as the optimized mapping in CROSS [5], decompose highprecision matrix operations into multiple low-precision substeps. The resulting partial outputs are then aggregated to recover the desired precision. However, this reconstruction/aggregation stage becomes a critical performance bottleneck, as it effectively serializes execution after the matrix multiplications. Moreover, because reconstruction is offloaded to vector processors outside the main matrix engine, it introduces data movement overhead and reduces systolic-array utilization. This inefficiency motivates our work: we integrate highprecision recomposition directly into the matrix multiplication dataflow. Our goal is to augment standard 8-bit systolic arrays to also perform the reconstruction step, producing full-precision matrix multiplication results without requiring computation outside the array. Crucially, we achieve this without relying on traditional reconfigurable SIMD-based multi-precision architectures [9], [14], which group adjacent processing elements to emulate wider arithmetic units. While effective in other domains, such spatial grouping incurs significant hardware overhead and reduces the effective array size available for parallel matrix multiplication as precision requirements increase. III. M ULTI -P RECISION M ATRIX M ULTIPLICATION IN E NHANCED S YSTOLIC A RRAYS TPU matrix engines (systolic arrays) are built around small, fixed-width multipliers and relatively wide accumulators. For instance, Google TPUs contain 8-bit integer multipliers in each processing element, while the column-wise accumulation (reduction) in a weight-stationary dataflow uses 32-bit adders [6]. This design makes them highly efficient for low-precision dense linear algebra. However, it also means that highprecision arithmetic, where input operands exceed 8 bits, cannot be executed directly within the array. Instead, larger integers must be decomposed into smaller digits, and matrix multiplication must be carried out over the corresponding low-precision elements. Intermediate results must then be aggregated using shift-and-add operations to reconstruct the final full-precision result. Systolic arrays are not designed to perform full-precision reconstruction internally; instead, this step is typically offloaded to vector processors outside the array that disrupts the otherwise uniform dataflow of matrix multiplication, introducing unnecessary data movement and effectively lowering the utilization of the matrix engine. To address this, we first show how full-precision input matrices can be decomposed into matrices composed of small digits (e.g., 8-bit integers). We then describe how standard matrix multiplication generates intermediate results, and finally how these results can be processed within the systolic array to recover the correct full-precision output as a seamless extension of the original dataflow.
12 47 35 68
×
54 23 71 89
=
3985 4459 6718 6857
b1 , b0 represent each digit. When we multiply them (a × b), you get a quadratic polynomial for: c = (a1 b1 ) · 102 + (a1 b0 + a0 b1 ) · 101 + (a0 b0 ) · 100
450320 2174 0 4 5 0 3 2 Decomposition × 5386 1 7 0 9 8 0 017098
2
We can express the calculation of these three coefficients as a matrix-vector multiplication a0 a1 0 b0 b1 = a0 b0 a1 b0 + a0 b1 a1 b1 0 a0 0
(low-precision matrix multiplication)
In this example and in the one shown in Fig. 1, the decomposed digits are intentionally placed least-significant digit first. This property is essential for proper operation scheduling within the systolic array, as explained in the next section. The multiplication of the single-digit matrices yields a 2 × 6 matrix of intermediate results (step 2 in Fig. 1). Recomposition then is required to combine the intermediate results to form the correct full-precision result. As shown in part 3 of Fig. 1 to perform this reconstruction step, neighbor columns corresponding to different digit positions should be combined through shifts and digit-level carry propagation thus completing the multi-digit addition. Carry propagation remains within the limits of each group of digits and does not expand to columns of different output digits.
=
1
15 67 33 28 99 57 69 99 34 87 137 54
Intermediate Result
3
15 69 67 99 +34 +33 4459 3985 28 87 99 137 +57 +54 6718 6857
Full-Precision Reconstruction
Fig. 1: Conversion of high precision matrix multiplication to low precision, through decomposition and reconstruction. A. Overview of the High-Precision Matrix Multiplication Flow Fig. 1 illustrates how a matrix multiplication with multidigit (high-precision) integers is transformed into an equivalent matrix multiplication over single-digit operands that can be computed efficiently on a systolic array. The example uses two-digit decimal numbers for clarity, but the same procedure applies to any other radix. The top of the figure shows the result of the multiplication of the full-precision matrices. Each full-precision matrix is first decomposed into an larger matrix consisting solely on smaller digits (step 1 in Fig. 1). For the left matrix, every element is expanded into a row vector containing the original digits but placed to consecutive columns from the least significant to most significant digit. This produces a wider matrix whose entries are single-digit values. The right matrix undergoes a different transformation. Instead of simply stacking digits, the constituent digits are arranged into a Toeplitz-like structure [15] that encodes all digit cross-products required by schoolbook digit-wise multiplication [16]. Each digit of the right-hand matrix is copied into multiple columns with appropriate horizontal shifts. These shifts ensure that when the single-digit matrices are multiplied, products between digit pairs automatically align with the correct positional weight. As a result, a single dense matrix multiplication generates all digit cross-products required by the original high-precision computation. To clarify this, let’s assume two entries of the corresponding input matrices, a is an element of left matrix and b an element of the write matrix. Since both contain two decimal digits we can write a = a1 ·10+a0 and b = b1 ·10+b0 , where a1 , a0 and
B. Proposed Systolic Array Architecture Performing precision reconstruction natively within systolic arrays without disrupting their typical dataflow requires almost no changes to their standard operation. We only need to augment the PEs of the final row of the systolic array with an additional adder, while the rest of the array remains unmodified. We base our architecture on a standard Google TPU-style systolic array, utilizing the widely adopted weight-stationary dataflow. In this baseline design, the right matrix (weights) is preloaded into the array, while the left matrix streams horizontally into the array from the western edge. Each PE performs an 8-bit multiplication between the incoming streaming data and its locally stored weight. The 16-bit products are then continuously accumulated vertically down each column using 32-bit adders. This typical structure is depicted in Fig. 2. The PEs of the last row called Reconstruction PEs (RPEs) include one additional input relative to the typical PE of the rest rows. These RPEs receive a digit from the RPE on the left and add it to the sum propagating downwards on the same column. In this way, digit carries are propagated horizontally along the final row of the array in sync with the vertical propagation of the per-column sums. The multiplexers stops this carry propagation by selecting a zero input to the adder when result exceeds the columns of the full-precision result. Figure 3 presents a cycle-by-cycle illustration of reconstruction performed concurrently with the systolic wavefront of matrix multiplication as it propagates through the array. Initially, the right-hand matrix obtained from the decomposition in Fig. 1 is assumed to be preloaded into the array, while the left-hand matrix is streamed in row by row from the left. To emphasize only on the operation of reconstruction, Fig. 3 considers only the first row of the decomposed left matrix from Fig. 1.
Fig. 2: The proposed systolic array. All PEs keep their typical multiply-add structure while the PEs of the last row (RPEs) are enhanced with an additional adder that combines horizontal digit-level carry propagation with vertical per-column addition. As the systolic array operates, the first row of the intermediate result matrix is progressively generated along the bottom of the array. The first outputs begin to appear in consecutive columns starting at cycle 3. Those elements effectively enter the RPEs of the last row. Each RPE retains the least significant byte for its column and forwards the most significant byte as a carry to the adjacent RPE on the right. This digit-level carry propagation continues until the full precision result has been reconstructed. The required number of digits is fixed at design time, with the appropriate configuration applied at runtime. Please note that the intermediate results are streamed in least-significant-digit-first order. This ordering is critical since it ensures that digit k arrives only after digit k−1 has produced its carry on the previous column. C. NTT-Specific Full-Precision Matrix Decomposition This section describes how full-precision matrix decomposition occurs specifically for NTT and simplified using the Basis-Aligned Transformation (BAT) [5]. The 4-step NTT factorizes an N -point transform into smaller transforms, computed as the product of an N × N input matrix A and a twiddle-factor matrix T . Suppose A and T contain b-bit full-precision integers. To execute this multiplication on the proposed systolic array, which only supports w-bit multipliers (w < b), both matrices are decomposed using the methodology from Section III-A. Each entry is decomposed into b k= w radix-2w digits. These digits correspond to positional weights 20 , 2w , . . . , 2(2k−2)w . This expansion structurally resizes matrix A from N × N to N × kN , and matrix T from N × N to kN × (2k − 1)N . Matrix T contains only twiddle factors that depend on the transform parameters (N, q) and are constant for a specific NTT instance. For this matrix, BAT suggests a way to reduce
Fig. 3: Running example on the in-sync execution of percolumn vertical reduction needed for single-digit matrix multiplication and horizontal digit-level carry propagation for recomposing full-precision matrix multiplication output. the size of this matrix taking into account that the contribution of each digit of decomposed matrix T should be weighted modulo q. Thus, all digit weights that exceed 2kw can be folded back to lower digits. Since matrix T is known at design time and constant at runtime this reallocation of digits can be performed offline. The following toy example illustrates the offline reduction
applied to twiddle factor digits. Consider the value 54 (the first element of the right matrix in Fig. 1) as one of the twiddle factors. Using base-10 digit decomposition with two-digit precision, the value is expanded into the following Toeplitz matrix: 0 1 2 10 10 10 4 5 0 54 → 0 4 5 Assuming that all computations are performed modulo 97 the contribution of digit 5 with weight 102 can be reduced offline to 5·102 mod 97 = 15 = 1·101 +5·100 . Distributing this result to the columns with weights 100 and 101 , respectively, yields 0 10 4 0
101 5 4
0 102 10 0 4 → 5 0+5
101 5 4 = 4+1 5
5 . 5
Twiddle-factor digit-folding is applied incrementally until all digits of the T matrix fit in w bits. After applying offline this redistribution of twiddle factors as suggested by BAT [5] the dimension of decomposed matrix T reduces to k N × k N . IV. E XPERIMENTAL E VALUATION To evaluate the proposed approach, we quantify both the performance speedup achieved when executing FHE NTT kernels (for transform sizes ranging from 212 to 216 ) using the 4-step NTT algorithm, and the minimal hardware overhead introduced relative to a standard TPU. For performance measurements, we utilize the cycle-accurate SCALE-Sim simulator [17], which models architectures similar to commercial AI accelerators. The simulated baseline features a 128 × 128 systolic matrix unit (MXU) paired with a 128-lane vector processing unit (VPU) [6]. Each MXU integrates dedicated 8 MB on-chip input, weight, and output buffers with sufficient bandwidth to sustain peak systolic throughput. These buffers deliver one operand per processing element per cycle (e.g., 128 8-bit words/cycle for a 128×128 array), ensuring continuous utilization of the array. We compare the proposed design with two state-of-theart approaches: (a) TensorFHE [3] executes low-precision matrix multiplications on the MXU and the full-precision reconstruction is offloaded to the VPU; (b) CROSS [5] follows the same approach as TensorFHE for matrix multiplication and result reconstruction but follows a different decomposition of the original matrices to their low-precision counterparts. The decomposition of CROSS is similar to the one described in Section III-A. In all cases, including the proposed approach, when matrix multiplication sizes exceeds the size of MXU, multiplication is performed in tiles [18]. Because FHE schemes dictate 32- to 64-bit modular arithmetic representing the most commonly used moduli in latticebased cryptography and typical FHE libraries, we evaluate both configurations. A. Runtime performance comparisons Table I summarizes the runtime comparison between TensorFHE, CROSS, and the proposed architecture for 32-bit and 64-bit full-precision inputs. For each NTT size, the table
TABLE I: Runtime comparison (in cycles) and relative speedup of the 4-step NTT algorithm between TensorFHE [3], CROSS [5] and the proposed architecture for 32-bit and 64-bit full-precision inputs across varying transform sizes. TensorFHE [3] CROSS [5] Proposed NTT Size MXU VPU Total MXU VPU Total MXU VPU Total cycles cycles cycles cycles cycles cycles 212 14240 2176 16416 3566 2176 5742 3566 256 3822 213 15264 4288 19552 9174 4288 13462 9174 448 9622 214 16288 8576 24864 16318 8576 24894 16318 896 17214 215 42816 17024 59840 42846 17024 59328 42846 1664 44510 216 81632 34048 115680 81662 34048 115710 81662 3328 84990
(a) 32-bit data width TensorFHE [3] NTT Size MXU VPU Total cycles cycles
CROSS [5] MXU VPU Total cycles cycles
Proposed MXU VPU Total cycles cycles
212 56960 8320 65280 14270 8320 22590 14270 256 14526 213 61056 16576 77632 36702 16576 53278 36702 448 37150 214 65152 33152 98304 65278 33152 98430 65278 896 66174 215 171264 66176 237440 171390 66176 237566 171390 1664 173054 216 326528 132352 458880 326654 132352 459006 326654 3328 329982
(b) 64-bit data width
reports the total number of cycles required to execute a complete 4-step NTT. For all architectures, the execution time is divided into MXU cycles, corresponding to 8-bit digit matrix multiplications, and VPU cycles, corresponding to all operations performed outside the systolic array. In TensorFHE and CROSS, VPU cycles include the shift-add operations required to reconstruct full-precision results, as well as modular reductions, twiddle-factor multiplications, and matrix transpositions involved in the 4-step NTT. In contrast, for the proposed architecture, VPU cycles account only for modular reduction, since full-precision reconstruction is performed online during matrix multiplication by the RPEs within the systolic array. Across both precision settings, the proposed architecture consistently reduces execution cycles, achieving speedups between 1.33× and 4.49×. The results for 32-bit and 64-bit operands exhibit nearly identical trends and speedups. The primary source of performance improvement is the elimination of VPU-based full-precision output reconstruction that involves shift-add operations that are not natively supported by typical MXUs. In TensorFHE [3] and CROSS [5], a non-negligible portion of execution time is spent outside the systolic array combining partial sums and performing digit-level carry propagation. By integrating these operations directly into the MXU dataflow, the proposed architecture removes this overhead entirely. The largest speedups are observed for smaller NTT sizes against TensorFHE [3]. For example, in the 32-bit configuration at size 212 , execution time is reduced from 16416 cycles to 3822 cycles, corresponding to a 4.30× improvement. In this regime, VPU execution constitutes a significant fraction of total runtime, so eliminating it yields substantial gains. As the NTT size increases, matrix multiplication increasingly dominates execution time. Consequently, the relative impact of removing VPU overhead decreases and speedups converge to 1.36× for the largest transform size 216 .
Fig. 4: The physical layout of an 8×8 instance of the proposed systolic array. TABLE II: Area and power comparison between a conventional systolic array and the proposed enhanced design with RPEs in the last row, implemented at 7 nm using the OpenROAD physical synthesis platform. Area (mm2 ) Power (W) baseline proposed overhead baseline proposed overhead 8×8 0.0041 0.0042 2.79% 0.0254 0.0241 5.52% 16 × 16 0.0165 0.0167 1.39% 0.0989 0.0963 2.76% 32 × 32 0.0659 0.0664 0.70% 0.3903 0.3850 1.38% 64 × 64 0.263 0.264 0.35% 1.84 1.86 0.69% 128 × 128 1.054 1.056 0.17% 7.39 7.41 0.34% 256 × 256 4.217 4.221 0.09% 29.56 29.61 0.17% MXU Size
B. Physical synthesis comparisons To evaluate the hardware overhead introduced by the extra adder and multiplexing logic in the RPEs, which replace the final row of standard PEs in the baseline systolic array, we implemented both the proposed design and the baseline using the OpenROAD physical synthesis platform with the ASAP7 technology library targeting a clock frequency of 1 GHz that matches industry-grade TPUs [6]. As an example, Fig. 4 shows the physical layout of a small-scale 8 × 8 instance of the proposed systolic array. Table II reports the post-layout area and power consumption for both architectures across multiple array sizes, alongside the relative overhead introduced by the additional full-precision output reconstruction logic. The results demonstrate that the proposed modifications incur minimal hardware penalties. Specifically, the area overhead decreases from 0.35% for small arrays to just 0.09% for larger configurations, while the power overhead drops from 0.69% to 0.17%. This scaling trend reflects the amortization of the final row’s additional logic as the overall array size increases. Overall, the proposed approach reduces the total energy per NTT by significantly lowering latency with only negligible power overhead, demonstrating that native full-precision output reconstruction within the systolic array is both an effective performance optimization and a practical energy efficient hardware enhancement.
V. C ONCLUSION Executing FHE kernels, such as the NTT, on AI ASICs offers a widely available, energy-efficient alternative to custom hardware. However, a fundamental discrepancy exists: AI ASICs are tuned for low-precision arithmetic, whereas FHE requires high precision. Mapping high-precision workloads onto low-precision matrix engines using standard dataflows is possible but creates serialization bottlenecks during fullprecision result reconstruction. To address this, we propose a novel multi-precision systolic array architecture that unifies low-precision matrix multiplication with in-situ high-precision reconstruction. Our design introduces a novel dataflow where matrix multiplication results flow vertically, perfectly synchronized with horizontal, digit-level carry propagation. This fused approach eliminates reconstruction bottlenecks, delivering significant execution speedups for NTT kernels with a hardware overhead of less than 1%. R EFERENCES [1] W. Wang et al., “Accelerating fully homomorphic encryption using gpu,” in IEEE Conf. on High Performance Extreme Computing, 2012, pp. 1–5. [2] N. Samardzic et al., “F1: A fast and programmable accelerator for fully homomorphic encryption,” in IEEE Inter. Symp. on Microarchitecture (MICRO), 2021, pp. 238–252. [3] S. Fan, Z. Wang, W. Xu, R. Hou, D. Meng, and M. Zhang, “TensorFHE: achieving practical computation on encrypted data using gpgpu,” in IEEE Inter. Symp. on High-Performance Comp. Arch. (HPCA), 2023, pp. 922–934. [4] Y. Sun et al., “Tensorfhe+: Fully homomorphic encryption acceleration based on linear algebra,” IEEE Trans. on Computers, vol. 75, no. 2, pp. 612–627, 2026. [5] J. Tong et al., “Leveraging ASIC AI chips for homomorphic encryption,” in IEEE Inter. Symp. on High Performance Comp. Arch. (HPCA), 2026. [6] N. P. Jouppi et al., “Ten lessons from three generations shaped google’s tpuv4i : Industrial product,” in IEEE Inter. Symp. on Comp. Arch. (ISCA), 2021, pp. 1–14. [7] A. Acar, H. Aksu, A. S. Uluagac, and M. Conti, “A survey on homomorphic encryption schemes: Theory and implementation,” ACM Computing Surveys, vol. 51, no. 4, pp. 1–35, 2018. [8] J.-C. Bajard, J. Eynard, M. A. Hasan, and V. Zucca, “A full RNS variant of fv like somewhat homomorphic encryption schemes,” in International Conf. on Selected Areas in Cryptography. Springer, 2016, pp. 423–442. [9] V. Camus, C. Enz, and M. Verhelst, “Survey of precision-scalable multiply-accumulate units for neural-network processing,” in IEEE Inter. Conf. on Artificial Intelligence Circ. and Syst. (AICAS), 2019, pp. 57–61. [10] T. Ajayi and D. Blaauw, “Openroad: Toward a self-driving, open-source digital layout implementation tool chain,” in Proceedings of Government Microcircuit Applications and Critical Technology Conf., 2019. [11] S. Halevi and V. Shoup, “Algorithms in helib,” in Annual Cryptology Conf. Springer, 2014, pp. 554–571. [12] D. Harvey and J. van der Hoeven, “Faster integer multiplication using short lattice vectors,” The Open Book Series, vol. 2, no. 1, pp. 293–310, 2019. [13] D. H. Bailey, “FFTs in external or hierarchical memory,” Journal of Supercomputing, vol. 4, no. 1, pp. 23–35, 1990. [14] W. Mao et al., “An energy-efficient mixed-bitwidth systolic accelerator for nas-optimized deep neural networks,” IEEE Trans. on Very Large Scale Integration (VLSI) Systems, vol. 30, no. 12, pp. 1878–1890, 2022. [15] R. M. Gray, “Toeplitz and circulant matrices: A review,” 2006. [16] D. E. Knuth, The art of computer programming: Seminumerical algorithms, volume 2. Addison-Wesley Professional, 2014. [17] R. Raj et al., “Scale-sim v3: a modular cycle-accurate systolic accelerator simulator for end-to-end system analysis,” in IEEE Inter. Symp. on Perf. Analysis of Systems and Software (ISPASS), 2025, pp. 186–200. [18] A. Schönhage and V. Strassen, “Schnelle multiplikation grosser zahlen,” Computing, vol. 7, no. 3, pp. 281–292, 1971.