A high-throughput multi-precision large integer multiplication optimization method and system based on GPU tensor core

Abstract

This invention discloses a high-throughput, multi-precision large integer multiplication optimization method and system based on GPU tensor cores, belonging to the field of cryptography. It addresses the issues of simultaneously guaranteeing computational accuracy and power in large integer multiplication in privacy computing applications, as well as low scalability and high cost. The method includes: targeting the variable-number × constant calculation pattern common in cryptographic batch processing, innovatively mapping large integer multiplication to a tensor core, supplemented by a continuity-aware rearrangement mechanism to eliminate carry bottlenecks, and a hardware alignment and packing mechanism to overcome memory access limitations. This securely reduces carry propagation to the registers and warps, and completely eliminates shared memory read bottlenecks, aiming to achieve perfect synergy between computational power and storage bandwidth without sacrificing computational accuracy. This invention is applicable to large-scale privacy computing applications such as fully homomorphic encryption (FHE) and zero-knowledge proofs (ZKP).

Description

Technical Field

[0001] This application relates to the field of cryptography, and in particular to high-throughput, multi-precision, large integer multiplication optimization based on GPU tensor cores. Background Technology

[0002] High-precision large integer multiplication is the most fundamental arithmetic primitive in privacy-preserving computing applications such as fully homomorphic encryption (FHE) and zero-knowledge proofs (ZKP). In these cryptographic systems, the core underlying operations are typically finite field arithmetic, with the main computational overhead coming from repeatedly performed large integer modular multiplication operations. As the parameter size increases, the bit width of the operands can reach thousands of bits, and the throughput of multi-precision multiplication directly determines the scalability of the end-to-end system.

[0003] In recent years, GPUs have been widely used to accelerate cryptographic workloads. However, current state-of-the-art GPU implementations of large number arithmetic (such as the CGBN library) primarily rely on traditional CUDA cores. This is because multi-precision multiplication has strict carry dependencies and irregular data flows. Meanwhile, to accelerate deep learning workloads that are primarily based on matrix multiplication (GEMM), modern GPUs have introduced Tensor Cores, which offer extremely high arithmetic throughput on low-precision integers (such as INT8 input, S32 accumulation). However, due to the mismatch between the carry dependencies and irregular data flows of large integer multiplication and the strict matrix multiplication semantics of Tensor Cores, existing GPU acceleration solutions have failed to effectively utilize the enormous computational power of Tensor Cores.

[0004] In practical cryptographic applications (such as Montgomery multiplication), there are often a large number of batch multiplication operations, where one operand (such as a modulus-related parameter) is a constant shared throughout the batch. This "variable × constant" computational model provides an opportunity to accelerate high throughput using tensor kernels.

[0005] Existing GPU multiprecision and large integer multiplication acceleration solutions primarily rely on traditional CUDA Cores (scalar arithmetic units), with representative implementations such as the widely used CGBN library. These solutions typically represent large integers using 32-bit or 64-bit limbs, performing multiplication through multiply-accumulate operations with carry-over. However, multiprecision integer multiplication inherently involves strict serialization carry dependencies and irregular data flows. Under traditional scalar architectures, this strong dependency severely limits instruction-level parallelism, making it difficult to fully utilize the GPU's peak computing power. As operand width increases (e.g., 4096 bits or 8192 bits), performance improvements in traditional solutions are often limited by complex instruction scheduling and cross-thread data exchange, making it difficult to achieve purely compute-bound efficient execution.

[0006] Another common acceleration approach is based on number-theoretic transforms (NTTs) or fast Fourier transforms (FFTs), which attempt to transform large integer multiplication into polynomial multiplication to reduce theoretical time complexity. While this method has theoretical advantages at very large bit widths, its complex memory access patterns often lead to severe memory-bound problems in medium-bit width scenarios (e.g., 2048 to 8192 bits) commonly used in cryptography and privacy computing. Furthermore, the high overhead of constant factors and the complexity of field parameter transformations make end-to-end integration of such schemes costly in practical low-level arithmetic systems, resulting in highly unstable throughput.

[0007] In recent years, with the evolution of hardware architecture, modern GPUs have introduced Tensor Cores specifically designed for matrix multiplication and accumulation (MMA), which provide orders of magnitude leaps in computational throughput for low-precision data types such as INT8. To take advantage of this hardware feature, some experimental methods map large integer multiplication to generalized matrix multiplication (GEMM) through byte decomposition and Toeplitz matrix transformation. However, directly adopting the standard GEMM mapping strategy faces two major obstacles: First, the carry chain breakage and misalignment problem. The standard high-performance matrix tiling strategy causes the computation results of a single warp to be fragmented in the logical output space. This forces carry operations with strong sequential dependencies to be performed frequently across warps, and the resulting high-latency shared memory communication will directly drag down the overall performance. Second, the severe memory access conflict problem. Tensor Core's dedicated data loading instructions (such as ldmatrix) are highly sensitive to 16-byte memory alignment, fixed mapping mode, and shared memory bank conflicts. The naive Toeplitz matrix construction not only wastes a lot of shared memory capacity, but also causes serious read conflicts, causing the tensor core pipeline to be in a waiting state.

[0008] Therefore, how to effectively resolve the inherent carry dependency and irregular memory access obstacles of large integer multiplication while fully unleashing the extremely high throughput of Tensor Cores has become the core challenge of current cryptographic underlying hardware acceleration design. Summary of the Invention

[0009] The purpose of this invention is to solve the technical problems of not being able to simultaneously guarantee the accuracy and computing power of large integer multiplication in privacy computing applications, as well as low scalability and high cost. It provides a high-throughput, multi-precision large integer multiplication optimization method and system based on GPU tensor cores.

[0010] The technical solution adopted by this invention to solve the above problems is: a high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores, the method comprising:

[0011] Step 1: After the kernel boots, it first creates operands for variables in the GPU shared memory. Allocate block cache space and invoke hardware-aligned shared memory packing mechanism to pack shared constant operands. The components are packaged into shared memory in reverse byte order, with redundant channel striping and staggered padding, for subsequent conflict-free reuse, and then thread bundle synchronization is performed.

[0012] Step 2: The calculation process is organized into nested loops; the outer loop iterates through the column-wise blocks of the output results, i.e., the logical weight dimension of large integers. The inner loop iterates and accumulates along the depth direction, i.e., the byte dimension after operand decomposition. ;

[0013] In the inner loop, the system continuously retrieves variables from global memory. The data blocks are loaded into shared memory and then into registers; simultaneously, the target column coordinates are calculated for each thread bundle and its channel according to the continuity-aware rearrangement strategy, resulting in the current Toeplitz shift at constant. The precise offset of the shared memory buffer; subsequently, the dedicated instruction ldmatrix is ​​invoked to load it conflict-free. The data fragment is written to the register, and the INT8 tensor core is called to execute the matrix multiply-add MMA instruction, and the result is accumulated in the S32 type register;

[0014] Step 3: When the inner layer After all the dimensional accumulations are completed, the intermediate S32 results stored in the registers of each thread are continuously distributed in the logical output dimension. The system directly calls the carry propagation algorithm of thread bundle parallelism, and converts the byte-level accumulation results into high-precision values ​​through register packing and shuffle instructions within the thread bundle, supplemented by shared memory boundary buffers, and writes the final results back to global memory.

[0015] Further, step 1 includes: large integer bit Represented as a polynomial with base 8:

[0016]

[0017]

[0018] and Multiplication can be written in convolution form:

[0019] .

[0020] Further, step 1 includes: constructing the Toeplitz structure matrix. :

[0021]

[0022] in, This represents the Toeplitz structure matrix constructed from the byte components of operand B, which is a structure matrix uniquely determined by B. for The element in the i-th row and t-th column;

[0023] Batch input Composition matrix , It is the set of all variable operands. Represents a set The number of elements, then the matrix multiplication result The result of a large integer multiplication before carry, the first... row representation The result before carry;

[0024]

[0025] right After performing carry processing on each row, the result of the large integer multiplication is obtained.

[0026] Furthermore, the hardware-aligned shared memory packing mechanism in step 1 specifically includes:

[0027] First, based on the byte length of the constant operand. The reference offset of the buffer required for a single channel is calculated in advance. and the span length of a single strip To provide sufficient contiguous addressing space for all convolutional shifts, the strip length is set to a constant length plus additional padding alignment space. ;

[0028] Subsequently, in order to naturally form the convolution shift effect of the Toeplitz matrix in physical memory, the system utilizes the multi-threaded concurrency capability of the GPU to reverse the original byte sequence of the constant operand B and preferentially write it to the base channel of shared memory, namely channel 0; on this basis, the system redundantly copies the reverse sequence to an additional 15 channel stripes.

[0029] During channel replication, the system applies a channel index to the starting address of each channel. Equal-length physical memory offsets; this alignment skew mechanism ensures that when the Tensor Core requests Toeplitz shift data for different columns according to the continuum rearrangement strategy, the actual physical address calculated for each channel is always a multiple of 16, thus meeting the 16-byte alignment check requirement of the ldmatrix.sync.aligned instruction.

[0030] Furthermore, the hardware-aligned shared memory packing mechanism in step 1 also includes: making additional offset adjustments to the storage locations of some channels to reduce the original 8-way bank conflicts to 4-way conflicts; the remaining conflicts are inherent to continuity-aware rearrangement.

[0031] Furthermore, the continuity-aware rearrangement strategy in step 2 includes:

[0032] First, define the matrix multiplication block size for large integer multiplication in the GPU shared memory, and set the thread bundle warp array for executing this block; for each thread, obtain its thread bundle index in the current computation block. and the channel index within the thread bundle ;

[0033] based on and Calculate the rearranged target column coordinates for each thread The specific calculation rules are implemented through cross-interleaving bitwise operations, and its mathematical mapping logic is as follows:

[0034]

[0035] in, To calculate the block size in the row direction within the block; This represents the block size in the column direction; within a block, the column indexes are rearranged by group. To indicate the size of the group;

[0036] According to the calculation Determine the load offset of the constant operand in shared memory. To meet the 16-byte alignment requirement of the underlying hardware loading instruction ldmatrix and to alleviate bank conflicts in shared memory, the system performs alignment and stepping processing on the offset.

[0037] Using the GPU's alignment matrix loading instructions, according to the calculated Constant operands are batch-loaded from shared memory into thread registers in rearranged order; then, the loaded variable operands and rearranged constant operands are sent into the tensor core to perform low-precision MMA atomic calculations, and the accumulated results are stored in high-precision registers.

[0038] Furthermore, the thread-beam parallel carry propagation algorithm in step 3 specifically includes:

[0039] First, pack the 8 bytes into two 32-bit values ​​locally;

[0040] The thread bundle synchronization instruction __shfl_up_sync is used to pass the carry generated by the last byte of the previous thread to the first byte of the next thread;

[0041] For carry operations that cross thread bundle boundaries, they are passed and accumulated through a shared memory buffer.

[0042] Secondly, the present invention provides a high-throughput multi-precision large integer multiplication optimization method system based on GPU tensor cores. The system has program modules corresponding to the steps of the high-throughput multi-precision large integer multiplication optimization method based on GPU tensor cores as described above, and executes the steps in the high-throughput multi-precision large integer multiplication optimization method based on GPU tensor cores during runtime.

[0043] Thirdly, the present invention provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, it executes the steps of the high-throughput multi-precision large integer multiplication optimization method based on GPU tensor cores as described above.

[0044] Fourthly, the present invention provides a computer-readable storage medium for storing a computer program that executes a high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores as described above.

[0045] The beneficial effects of this invention are:

[0046] This invention provides a high-throughput, multi-precision integer multiplication method and system specifically for GPU Tensor Cores. Addressing the "variable × constant" (e.g., Montgomery reduction) computation pattern common in cryptographic batch processing, it employs an innovative continuity-aware permutation and redundant packing memory management mechanism to securely reduce carry propagation to within registers and warps, completely eliminating shared memory read bottlenecks. This invention aims to achieve perfect synergy between computing power and storage bandwidth without sacrificing computational accuracy, providing a robust and highly scalable arithmetic foundation for large-scale privacy computing applications such as fully homomorphic encryption (FHE) and zero-knowledge proofs (ZKP).

[0047] This invention is applicable to large-scale privacy computing applications such as fully homomorphic encryption (FHE) and zero-knowledge proof (ZKP), and is particularly suitable for batch processing workloads with shared operands. Attached Figure Description

[0048] To more clearly illustrate the technical solution of this application, the drawings used in the embodiments will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0049] Figure 1 Here is a flowchart of the byte decomposition and transformation process;

[0050] Figure 2 MMA atoms;

[0051] Figure 3 This is a combination of traditional MMA solutions;

[0052] Figure 4 For continuity-aware rearrangement strategies;

[0053] Figure 5 Loading process for ldmatrix;

[0054] Figure 6 For redundant packing layout of operand B in shared memory;

[0055] Figure 7 The results are from the throughput-batch size experiment.

[0056] Figure 8 The results are experimental results with relative throughput peak normalization. Detailed Implementation

[0057] The specific implementation method described herein, a high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores, includes:

[0058] The goal of this implementation is to design a high-throughput, multi-precision multiplication primitive specifically for GPU tensor cores, particularly suitable for batch processing workloads with shared operands. This system reconstructs byte-level multiplication to match tensor core semantics while maintaining carry locality and satisfying memory access constraints. The overall architecture is based on the following four core mechanisms:

[0059] 2.1 Byte Decomposition and Toeplitz Matrix Transformation

[0060] To leverage the high-throughput INT8 tensor cores on modern GPUs, this implementation decomposes large integer operands into polynomial representations with a base of 2^8. The product of two large integers can be equivalent to a discrete convolution over a sequence of bytes: large integer bit It can be represented as a polynomial with a base of 8.

[0061]

[0062]

[0063] So, and Multiplication can be written in the form of convolution as follows:

[0064]

[0065] To map this convolutional structure to Generalized Matrix Multiplication (GEMM), this implementation constructs a Toeplitz matrix using shared operands. For batch operations common in cryptographic workloads, multiple independent variable operands are aggregated into a dense matrix, thereby converting irregular convolutions into standard U8×U8→S32 matrix multiplication. This transformation allows computation to be performed efficiently directly on the tensor kernel.

[0066] Specifically, constructing the Toeplitz structure matrix :

[0067]

[0068] Batch input Composition matrix The result of matrix multiplication is... The result of a large integer multiplication before carry, the first... row representation The result before carry, batch calculation of multiple different large integers. Multiply by the same large integer B, Let A represent the set of all distinct large integers A.

[0069]

[0070] right After performing carry-over processing on each row, the result of the large integer multiplication is obtained. The overall transformation process is as follows: Figure 1 As shown.

[0071] 2.2 Continuity-Aware Rearrangement Strategy

[0072] While the Toeplitz transformation formally conforms to the GEMM definition, the standard GEMM block strategy is not suitable for large integer multiplication. In the standard Tensor Core execution model, to maximize arithmetic strength, multiple matrix multiply-add (MMA) atomic operations are typically aggregated into larger blocks and distributed across multiple warps. This results in the computational results produced by a single warp being fragmented in the logical output space, forcing carry propagation to traverse different warps and relying on slow shared memory for communication.

[0073] Specifically, on NVIDIA graphics card Tensor Cores, matrix multiplication is performed in blocks. The smallest block that completes a single matrix multiplication is an execution unit called a thread bundle. We call a single thread bundle performing one matrix multiplication an MMA atom. Figure 2 This demonstrates the distribution of registers of various threads within an MMA atom in a matrix.

[0074] In real-world scenarios, to improve the computation / read ratio and fully utilize the computing power of Tensor Cores, conventional solutions often combine multiple thread bundles and perform multiple MMA atomic calculations within a single matrix multiplication block. For example... Figure 3 As shown, the two thread bundles each performed two MMA atomic calculations. It can be observed that the calculation results of each thread bundle are separated by the other thread bundle, which means that subsequent carry propagation can only use shared memory to exchange data, rather than using an efficient intra-thread data exchange scheme (shuffle).

[0075] To reconcile the parallelism of the tensor core with the serial dependency of the carry chain, this implementation introduces a continuity-aware permutation strategy, such as... Figure 4 As shown. This strategy makes the results of multiple MMA atomic operations within a single thread bundle logically continuous by reordering the memory access patterns of shared operands. This continuity allows the carry chain to be resolved primarily through local register operations and efficient in-thread shuffle instructions, significantly reducing reliance on shared memory.

[0076] In order to achieve the rearrangement, the system in this embodiment performs the following steps:

[0077] First, define the matrix multiplication block size for large integer multiplication (e.g., BLK_M, BLK_N, BLK_K) in the GPU shared memory, and set the warp array (e.g., LAYOUT_M × LAYOUT_N) for executing this block. For each thread, obtain its warp index in the current computation block (…). ) and the channel index within the thread bundle ( (Value range: 0~31).

[0078] To ensure the results of multiple matrix multiply-add (MMA) atomic operations within a single thread bundle are logically consecutive, breaking the fragmented layout of native MMA, the system is based on... and Calculate the rearranged target column coordinates for each thread The specific calculation rules are implemented through cross-interleaving bitwise operations, and its mathematical mapping logic is as follows:

[0079]

[0080] in, , , Preset system configuration parameters for controlling the rearrangement granularity and thread bundle arrangement.

[0081] The system calculates based on the previous step Determine the load offset of the constant operand in shared memory. To meet the 16-byte alignment requirement of the underlying hardware loading instruction ldmatrix and to alleviate bank conflicts in shared memory, the system performs offset alignment and additional stepping.

[0082] Next, the system uses GPU-specific alignment matrix loading instructions to load the matrix according to the calculated... Constant operands are batch-loaded from shared memory into thread registers in rearranged order. Subsequently, the loaded variable operands and the rearranged constant operands are sent to the tensor core to perform low-precision (e.g., INT8) MMA atomic calculations, and the accumulated results are stored in high-precision (e.g., S32) registers.

[0083] After the aforementioned reordering loading and computation, the S32 accumulated result blocks output by the Tensor Core to each thread's register automatically exhibit strict physical continuity in the logical output columns (i.e., the consecutive byte weights of large integers). This eliminates the need to write intermediate results back to shared memory for cross-thread bundle reassembly during subsequent large integer carry propagation.

[0084] 2.3 Hardware-aligned shared memory packing mechanism

[0085] In the actual loading of the Toeplitz matrix, if it is explicitly generated in GPU shared memory... Large Toeplitz matrices not only waste valuable shared memory but also lead to severe memory access conflicts. Modern Tensor Core's high-performance pipelines heavily rely on specific aligned load instructions (such as ldmatrix.sync.aligned.m8n8) to efficiently move shared memory data into registers. These instructions have two stringent hardware constraints: first, the base address for each load must strictly adhere to 16-byte alignment; second, the instructions consistently allocate specific thread lanes to retrieve specific matrix columns, and improper layout can lead to severe read serialization (bank conflict) if multiple lanes simultaneously access the same memory bank. Figure 5 This demonstrates the underlying mapping characteristics of the Tensor Core Loading Instruction (ldmatrix). For operand B, each participating LoaderLane is actually responsible for retrieving an entire column of data from the Tensor Multiplication Atom (MMA Atom). This fixed mapping is the hardware basis for the shared memory packing strategy designed in this implementation.

[0086] To efficiently supply data while meeting the aforementioned hardware constraints, this implementation devises a lane-aware redundant packing mechanism, such as... Figure 6 As shown, the system does not store dense Toeplitz matrices. Instead, it allocates a precisely computed, compact buffer that transforms constant operands (large integers B) into multiple "lane stripes" with specific offsets.

[0087] Specifically, the system first determines the length of the constant operand in bytes. The reference offset of the buffer required for a single channel is calculated in advance. ) and the span length of a single strip ( To provide sufficient contiguous addressing space for all convolutional shifts, the strip length is set to a constant length plus additional padding alignment space. ).

[0088] Subsequently, in order to naturally form the convolution shift effect of the Toeplitz matrix in physical memory, the system utilizes the multi-threaded concurrency capabilities of the GPU to reverse the original byte sequence of the constant operand B and preferentially write it to the base channel (i.e., channel 0) of shared memory. Based on this, the system redundantly copies this reverse sequence to an additional 15 channel stripes.

[0089] During channel replication, the system applies a channel index to the starting address of each channel. Equal-length physical memory offsets. This crucial alignment skew mechanism ensures that when the tensor core requests Toeplitz shift data from different columns according to a contiguous rearrangement strategy, the actual physical address calculated for each channel is always a multiple of 16, thus perfectly matching the 16-byte alignment check requirement of the ldmatrix.sync.aligned instruction.

[0090] Furthermore, to mitigate bank conflicts, this implementation makes additional offset adjustments to the storage locations of some channels (e.g., an offset of 32 bytes), reducing the original 8-way bank conflicts to 4-way conflicts. The remaining conflicts are inherent to continuity-aware rearrangement, but this level of mitigation is sufficient to restore most of the memory bandwidth while preserving a carry-friendly output layout.

[0091] 2.4 Carry Propagation Algorithm with Thread Bundle Parallelism

[0092] The GEMM phase generates intermediate sums distributed across the various thread bundle register files. To compute the final result, carry propagation must be performed to convert these radix 2^8 fragments into standard radix 2^32 values ​​(limbs).

[0093] Thanks to the continuity-aware reordering strategy, each thread holds a logically continuous stripe of results (e.g., 8 consecutive byte accumulators). This implementation achieves an efficient carry chain:

[0094] 1. Vector Packing: First, pack the 8 bytes into two 32-bit values ​​locally.

[0095] 2. Intra-thread propagation: Using the efficient thread bundle synchronization instruction (__shfl_up_sync), the carry generated by the last byte of the previous thread is passed to the first byte of the next thread.

[0096] 3. Propagation between thread bundles: For carry that crosses the boundary of a thread bundle, it is passed and accumulated through a shared memory buffer.

[0097] Algorithm 1 demonstrates the complete carry propagation algorithm.

[0098]

[0099]

[0100] 2.5 Overall System Execution Flow

[0101] The large integer multiplication optimization method proposed in this implementation deeply integrates the above-mentioned independent core mechanisms (byte decomposition and Toeplitz transformation, continuity-aware rearrangement, hardware-aligned shared memory packing, and thread bundle parallel carry) to build a unified high-performance GPU Tensor Core kernel execution pipeline.

[0102] The execution of the entire algorithm kernel is mainly divided into three stages:

[0103] Preprocessing and initialization phase: After the kernel starts, it first prepares the variable operands in the GPU shared memory. Allocate tiled cache space and invoke hardware-aligned shared memory packing mechanism to pack shared constant operands. The data is packaged into shared memory in one go using reverse byte order, redundant channel striping, and staggered padding for subsequent conflict-free reuse. Thread bundle synchronization is then performed.

[0104] Block Tensor Multiply-Accumulate (MMA) Phase: The computation process is organized as nested loops. The outer loop iterates through the column-wise blocks of the output results (i.e., the logical weight dimension of large integers). The inner loop iterates through and accumulates the depth direction (i.e., the byte dimension after operand decomposition). The cumulative depth direction refers to the shared reduction dimension of GEMM, that is, the K-dimensional dimension common to A_mat (M x K) and T(B)(K x 2K). Therefore, the inner loop actually traverses the sub-blocks along this K-dimensional dimension. Within the inner loop, the system continuously loads variables from global memory. The data blocks are loaded into shared memory and then into registers; simultaneously, the target column coordinates are calculated for each thread bundle and its channel according to the continuity-aware rearrangement strategy, resulting in the current Toeplitz shift at constant. The precise offset of the shared memory buffer. Then, a dedicated instruction (ldmatrix) is invoked to load it conflict-free. The data fragment is written to the register, and the INT8 tensor core is called to execute the matrix multiply-add (MMA) instruction, which accumulates the result in an S32 type register.

[0105] Carry-and-merge / write-back phase: When the inner layer After all the dimensional accumulations are completed, the intermediate results stored in the S32 registers (i.e., the 32-bit signed integer registers of each thread) exhibit a strictly continuous distribution along the logical output dimension. The system directly calls the carry propagation algorithm in thread-bundle parallelism, using register packing and thread-bundle shuffle instructions, supplemented by a small amount of shared memory boundary buffers, to quickly convert the byte-level accumulation results into standard high-precision values, and safely write the final results back to global memory.

[0106] The formal definition of the complete execution pipeline is shown in Algorithm 2.

[0107]

[0108]

[0109] 3. Effects of the Invention

[0110] 3.1 Experimental Platform

[0111] This invention was validated on two different modern NVIDIA GPU architectures to demonstrate its universality in consumer and data center hardware. Table 1 shows the specific hardware specifications and environment configuration of the validation platform.

[0112] Table 1 Specifications of the Verification Platform

[0113]

[0114] This invention selects a typical batch processing scenario of large integer multiplication in cryptographic applications (i.e. ,in The batch size (BatchSize) is evaluated for shared constants and is set to a default value. The evaluation covered different operand widths from 2048 bits, 4096 bits, 6144 bits to 8192 bits.

[0115] To verify the advancement of this invention, the following three mainstream large integer arithmetic baseline systems were selected for comparison:

[0116] 1. Base: Standard baseline CUDA implementation, without using Tensor Cores for regular multiprecision multiplication.

[0117] 2. CGBN: NVIDIA's official open-source Cooperative Groups Big Number library, is currently the most widely used standard implementation of large number operations for high-performance GPU scalar cores (CUDA Cores) in both industry and academia.

[0118] 3. NTT: The state-of-the-art open-source GPU implementation of large number arithmetic based on number theoretic transformation, representing an alternative to processing extremely large integers by optimizing algorithmic complexity.

[0119] This invention selects execution time and throughput (MOps / s) as basic performance metrics, and uses Nsight Compute to analyze tensor core utilization and microarchitectural pauses as deeper hardware metrics. Simultaneously, its end-to-end throughput in Montgomery multiplication was tested.

[0120] 3.2 Experimental Results

[0121] The present invention conducted multi-dimensional performance evaluations under the above-mentioned verification platform and baseline settings. The experimental results mainly cover basic throughput, batch processing scalability, microarchitecture execution efficiency, ablation experiments, and end-to-end performance.

[0122] 3.2.1 Throughput Performance

[0123] Batch size is Under the given settings, the performance of the present invention and the baseline system at different bit widths is shown in the figure.

[0124] Table 2 Execution Time (milliseconds)

[0125]

[0126] The experimental results are shown in Table 2. As the operand width increases, the workload gradually shifts from being limited by memory bandwidth to being limited by computational power. The tensor core acceleration architecture designed in this invention demonstrates an overwhelming advantage. In 8192-bit large integer multiplication, the execution time of this invention is only 13.03 ms on an RTX 4080 Super and only 11.24 ms on an A800.

[0127] Table 3 Throughput (MOps / s)

[0128]

[0129] As shown in Table 3, the throughput metrics of this invention reached 80.47 MOps / s on the RTX 4080 Super and 93.29 MOps / s on the A800. Compared to the current state-of-the-art CGBN library (35.48 MOps / s on the A800), this invention achieved significant speedups of 1.75x and 2.63x, respectively. Furthermore, compared to the NTT-based implementation (only 7.91 MOps / s on the A800), this invention achieved a performance improvement of over 10x. This fully demonstrates the success of this invention in... The complexity of large number multiplication is efficiently mapped to the high-performance pipeline of the INT8 tensor core.

[0130] To evaluate the system's resource utilization efficiency under different loads, this invention fixes the operand width at 8192 bits and reduces the batch size from... Expand to The experimental results are as follows: Figure 7 and Figure 8 As shown.

[0131] Experiments show that the method of this invention not only has extremely high peak throughput but also excellent rapid saturation characteristics. This is particularly evident in relatively small batch sizes (e.g., ...). At this scale, the present invention achieved a throughput of 25.60 MOps / s, which is nearly 2.9 times that of CGBN (8.77 MOps / s) at the same scale. More importantly, the present invention achieves a throughput of 25.60 MOps / s at a batch size of 25.60 MOps / s. The system already reached 94.6% of its peak performance, while the CGBN and baseline systems only achieved about 60%-70% of their performance. This demonstrates that the byte decomposition and thread bundle parallelism strategy of this invention can effectively hide instruction and memory latency, significantly reducing the dependence on a large number of concurrent threads.

[0132] Microarchitecture profiling analysis using NVIDIA NsightCompute further reveals the fundamental source of the performance advantage:

[0133] 1. Tensor core utilization leaps: The TensorPipe utilization of this invention is as high as 62.70%, while all baseline methods (Base, CGBN, NTT) have a metric of 0% because they rely entirely on traditional CUDA cores.

[0134] 2. Significantly Reduced Memory and Scheduling Pauses: Traditional CGBN solutions exhibit high dependency wait pauses (StallWait reaches 32.84%), while NTT solutions are deeply mired in memory bottlenecks (DRAM utilization 59.03%, StallScoreboard 43.86%). In contrast, the pauses in this invention are mainly concentrated in the mathematical computation phase (StallMath accounts for 27.84%), meaning that the computation pipeline is executing arithmetic instructions at full load most of the time. This demonstrates that the "continuity-aware rearrangement" and "hardware-aligned packing" mechanisms in this invention successfully transform traditional memory / control-bound workloads into ideal compute-bound workloads.

[0135] 3.2.2 Montgomery multiplication end-to-end performance

[0136] This invention not only provides high-throughput basic multiplication primitives, but also directly applies them to the core component of finite field arithmetic—Montgomery multiplication (REDC algorithm)—and performs a complete end-to-end performance evaluation. Experiments were conducted to test... The experimental results of the independent 8192-bit Montgomery multiplications are shown in Table 4.

[0137] Table 4 Montgemoery multiplication throughput (KOps / s)

[0138]

[0139] After using the TC-Mult of this invention as the underlying acceleration engine, the system throughput reached 19,405 KOps / s. In comparison, the native REDC implementation of the CGBN library only achieved 4,614 KOps / s (this invention achieves approximately 4.2 times the native algorithm-level speedup); even compared to CGBN's highly optimized CIOS (coarse-grained integrated operand scan) implementation (12,419 KOps / s), this invention still achieves an absolute end-to-end speedup of 1.56 times.

[0140] Experimental results show that this invention achieves high tensor core pipeline utilization, successfully transforming this traditionally data-dependent primitive into a computationally intensive task on modern GPUs.

[0141] 3.2.3 Experiment Summary

[0142] Based on the comprehensive evaluation results from multiple dimensions, it can be concluded that the present invention significantly outperforms mainstream solutions represented by CGBN in terms of basic throughput, resource utilization efficiency, and actual end-to-end application performance.

[0143] By innovatively mapping large integer multiplication to a tensor core GEMM, supplemented by a continuity-aware rearrangement mechanism to eliminate carry bottlenecks, and a hardware alignment and packing mechanism to overcome memory access limitations, this invention successfully overcomes the hardware compatibility barrier that has long constrained GPU multi-precision computing. With up to 2.6x speedup for multiplication primitives and over 1.56x speedup for cryptographic end-to-end computation, this invention fully demonstrates its superior performance and stability in high-concurrency, large-bit-width computing scenarios, providing a practical computing foundation for building next-generation high-performance secure computing systems such as zero-knowledge proofs and fully homomorphic encryption.

Claims

1. A high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores, characterized in that, The method includes: Step 1: After the kernel boots, it first creates operands for variables in the GPU shared memory. Allocate block cache space and invoke hardware-aligned shared memory packing mechanism to pack shared constant operands. The components are packaged into shared memory in reverse byte order, with redundant channel striping and staggered padding, for subsequent conflict-free reuse, and then thread bundle synchronization is performed. Step 2: The calculation process is organized into nested loops; the outer loop iterates through the column-wise blocks of the output results, i.e., the logical weight dimension of large integers. The inner loop iterates and accumulates along the depth direction, i.e., the byte dimension after operand decomposition. ; In the inner loop, the system continuously retrieves variables from global memory. The data blocks are loaded into shared memory and then into registers; simultaneously, the target column coordinates are calculated for each thread bundle and its channel according to the continuity-aware rearrangement strategy, resulting in the current Toeplitz shift at constant. The precise offset of the shared memory buffer; subsequently, the dedicated instruction ldmatrix is ​​invoked to load it conflict-free. The data fragment is written to the register, and the INT8 tensor core is called to execute the matrix multiply-add MMA instruction, and the result is accumulated in the S32 type register; Step 3: When the inner layer After all the dimensional accumulations are completed, the intermediate S32 results stored in the registers of each thread are continuously distributed in the logical output dimension. The system directly calls the carry propagation algorithm of thread bundle parallelism, and converts the byte-level accumulation results into high-precision values ​​through register packing and shuffle instructions within the thread bundle, supplemented by shared memory boundary buffers, and writes the final results back to global memory.

2. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 1, characterized in that: Step 1 includes: large integer bit Represented as a polynomial with base 8: and Multiplication can be written in convolution form: 。 3. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 2, characterized in that: Step 1 includes: constructing the Toeplitz structure matrix. : in, This represents the Toeplitz structure matrix constructed from the byte components of operand B, which is a structure matrix uniquely determined by B. for The element in the i-th row and t-th column; Batch input Composition matrix , It is the set of all variable operands. Represents a set The number of elements, then the matrix multiplication result The result of a large integer multiplication before carry, the first... row representation The result before carry; right After performing carry processing on each row, the result of the large integer multiplication is obtained.

4. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 1, characterized in that: The hardware-aligned shared memory packing mechanism in step 1 specifically includes: First, based on the byte length of the constant operand. The reference offset of the buffer required for a single channel is calculated in advance. and the span length of a single strip To provide sufficient contiguous addressing space for all convolutional shifts, the strip length is set to a constant length plus additional padding alignment space. ; Subsequently, in order to naturally form the convolution shift effect of the Toeplitz matrix in physical memory, the system utilizes the multi-threaded concurrency capability of the GPU to reverse the original byte sequence of the constant operand B and preferentially write it to the base channel of shared memory, namely channel 0; on this basis, the system redundantly copies the reverse sequence to an additional 15 channel stripes. During channel replication, the system applies a channel index to the starting address of each channel. Equal-length physical memory offsets; this alignment skew mechanism ensures that when the Tensor Core requests Toeplitz shift data for different columns according to the continuum rearrangement strategy, the actual physical address calculated for each channel is always a multiple of 16, thus meeting the 16-byte alignment check requirement of the ldmatrix.sync.aligned instruction.

5. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 1, characterized in that: The hardware-aligned shared memory packing mechanism in step 1 also includes: making additional offset adjustments to the storage locations of some channels to reduce the original 8-way bank conflicts to 4-way conflicts; the remaining conflicts are inherent to continuity-aware rearrangement.

6. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 1, characterized in that: The continuity-aware rearrangement strategy in step 2 includes: First, define the matrix multiplication block size for large integer multiplication in the GPU shared memory, and set the thread bundle warp array for executing this block; for each thread, obtain its thread bundle index in the current computation block. and the channel index within the thread bundle ; based on and Calculate the rearranged target column coordinates for each thread The specific calculation rules are implemented through cross-interleaving bitwise operations, and its mathematical mapping logic is as follows: in, To calculate the block size in the row direction within the block; This represents the block size in the column direction; within a block, the column indexes are rearranged by group. To indicate the size of the group; According to the calculation Determine the load offset of the constant operand in shared memory. To meet the 16-byte alignment requirement of the underlying hardware loading instruction ldmatrix and to alleviate bank conflicts in shared memory, the system performs alignment and stepping processing on the offset. Using the GPU's alignment matrix loading instructions, according to the calculated Constant operands are batch-loaded from shared memory into thread registers in rearranged order; then, the loaded variable operands and rearranged constant operands are sent into the tensor core to perform low-precision MMA atomic calculations, and the accumulated results are stored in high-precision registers.

7. The high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores according to claim 1, characterized in that: The thread-beam parallel carry propagation algorithm in step 3 specifically includes: First, pack the 8 bytes into two 32-bit values ​​locally; The thread bundle synchronization instruction __shfl_up_sync is used to pass the carry generated by the last byte of the previous thread to the first byte of the next thread; For carry operations that cross thread bundle boundaries, they are passed and accumulated through a shared memory buffer.

8. A high-throughput, multi-precision, large integer multiplication optimization method system based on GPU tensor cores, characterized in that: The system has a program module corresponding to the steps of any one of claims 1-7, and executes the steps in the high-throughput multi-precision large integer multiplication optimization method based on GPU tensor cores when it runs.

9. A computer device, characterized in that: It includes a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, it performs the steps of the high-throughput multi-precision large integer multiplication optimization method based on GPU tensor cores as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that: The storage medium is used to store a computer program that executes the high-throughput, multi-precision, large integer multiplication optimization method based on GPU tensor cores as described in any one of claims 1-7.

Images