Conceptio › Archive › arXiv CS

arXiv CSopen access

Leveraging SIMD for Accelerating Large-number Arithmetic

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

Leveraging SIMD for Accelerating Large-number Arithmetic Subhrajit Das

Indian Institute of Technology, Gandhinagar India [email protected]

Abhishek Bichhawat

Indian Institute of Technology, Gandhinagar India [email protected]

5× more time than using a 2048-bit key [17]. While Intel and AMD are progressively investing more die area to support wider SIMD units to expand the horizon of application utilizing the data parallelism [1, 3], workloads involving largenumber arithmetic operations are falling behind in terms of performance due to their reliance on the scalar pipeline. Large-number arithmetic operations carry strong sequential dependencies between data elements, or limbs. For instance, when adding two 256-bit numbers (or 4 limbs on a 64-bit machine), each 64-bit limb is added separately, and the generated carry-bit is propagated to the next limb. As SIMD does not propagate the carry-bit across data elements in the same run, the overhead increases. Prior work [37, 69, 82] attempted to break these dependencies, but with limited success. As we show later, when adding two large numbers, managing carry-bits still requires 9× to 12× more time than the actual addition operation [69, 82]. Similarly, multiplying two large numbers [37] introduces long read-after-write chains, resulting in substantial throughput reduction. Through quantitative and qualitative analysis, we observe that the state-of-the-art approaches directly apply the existing sequential algorithm to SIMD, thereby inheriting the original dependency structure. Furthermore, we observe that the SIMD computation itself is never the performance bottleneck. Instead, the overhead of preparing and routing data for the intermediate SIMD operations imposed by the algorithm’s dependencies erodes the potential benefits of SIMD. In this paper, our objective is to address these shortcomings by asking the question: Can we reconstruct the algorithms to compute the more independent, data-parallel operations separately from the dependent operations, thereby improving effective SIMD utilization? We answer this question by proposing DigitsOnTurbo (DoT), an approach that leverages SIMD for accelerating large-number arithmetic operations. As the major bottleneck for addition and subtraction is the propagation of carry-bits across limbs, we divide the operation into four phases, such that the more common but faster tasks like limb addition can be done in parallel, while isolating the more performance-heavy but rarer task of handling cascading carries. For multiplication, we use the vertical and crosswise multiplication technique [57, 62, 80], which allows us to perform independent crosswise multiplication operations before dependencies are accounted for.

arXiv:2604.21566v1 [cs.DC] 23 Apr 2026

Abstract Large-number arithmetic, widely used in scientific computing and cryptography, has seen limited adoption of single instruction, multiple data (SIMD) parallelism on modern CPUs due to the inherent dependencies in traditional algorithms. We present DigitsOnTurbo (DoT), which restructures the computation around independent, data-parallel operations, rather than vectorizing the standard algorithms, thereby leveraging the benefits provided by SIMD. Over prior SIMD implementations, DoT achieves up to 1.85× speedups for addition and subtraction, and 2.3× for multiplication. When integrated into state-of-the-art libraries, DoT yields up to 4× speedup for addition and subtraction, and up to 2× speedup for multiplication, cascading into end-to-end throughput gains of up to 19.3% for scientific computations, and up to 7.9% latency and 5.9% throughput improvements on cryptographic implementations.

Yuvraj Patel

University of Edinburgh United Kingdom [email protected]

Introduction

Single instruction, multiple data (SIMD) allows modern CPUs to exploit data-parallel execution by running a single (same) instruction on multiple data elements, simultaneously [29, 39]. Different applications like linear algebra libraries (e.g., OpenBLAS [64], Intel MKL [45]), and computer vision and AI/ML frameworks [2, 5] have adopted SIMD by using specialized instructions to maximize CPU throughput. Interestingly, this has not been the case with applications performing large-number arithmetic that involve operations on operands that range from hundreds to thousands of bits. Large-number arithmetic is extensively used in high-precision scientific computing [7–9, 41, 60, 71, 81], and cryptographic applications [10, 22, 27, 31, 38, 40, 49, 61, 63, 66, 70, 75]; however, the libraries used for performing these operations do not capitalize on the benefits provided by SIMD architectures. The GNU multiple-precision arithmetic library (GMP) [34], one of the most widely deployed arbitrary-precision libraries, avoids the use of SIMD for accelerating these operations because “[SIMD does not provide] much support for propagating the sort of carries that arise in GMP” [73]. The cost of operating on large numbers shows up directly in workloads that depend on these primitives. For instance, computing 𝜋 to 1𝑀 digits takes 18× longer than for 100𝐾 digits [33], and decryption using a 4096-bit RSA key takes 1

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

We implement DoT in C with x86-64 SIMD intrinsics [47]. For addition and subtraction, DoT reduces the carry management overhead by over 2× compared with prior SIMD approaches [69, 82], delivering effective speedups of up to 1.85× over both those works and an optimized scalar addwith-carry baseline on an Intel Emerald Rapids CPU. DoT multiplication achieves a 2.31× speedup over prior work [37]. We integrate DoT into the GMP library (DoTMP) and OpenSSL (DoTSSL) by replacing the libraries’ addition and subtraction primitives and the base case multiplication operation. On end-to-end benchmarks, DoTMP improves GMPbench’s [33] overall score by 7.8%, with gains cascading into other operations: division improves by 8.4%, 𝜋 computation by up to 19.3%, and multiplication by up to 48.7% for the larger operands. DoTSSL improves OpenSSL throughput by up to 5.9% for FFDH and 5.2% for DSA, and reduces latency by up to 7.9 % for DSA, 6.1 % for FFDH, and 5.5 % for RSA. In summary, the contributions of this work are: • DoT add/sub: A 4-phase algorithm that reduces carry management overhead by nearly 2× and delivers strong SIMD speedups (1.8×), with the carry adjustment phase provably negligible for random inputs. • DoT multiplication: A SIMD multiplication routine that exposes all partial products as independent operations, eliminating RAW hazards and delivering 2× speedups for 256-bit multiplication. • End-to-end impact: DoTMP improves GMPbench’s overall score by 7.8% (𝜋 up to +19.3%, multiplication up to +48.7%); DoTSSL improves RSA, FFDH, and DSA throughput by up to 5.9% and reduces latency by up to 7.9%. All gains come from 1, 013 lines of portable C with intrinsics, ∼40 lines of integration changes across GMP and OpenSSL, and no hand-written assembly.

Background and Motivation

2.1

Large numbers and their representation

to maximize bits per limb [27, 69], while an unsaturated (reduced-radix) representation keeps headroom in each limb (e.g., base 251 for P-256 [36]) to simplify carry handling and field-specific reduction. 2.2

SIMD and Large-number Arithmetic

On x86-64, SIMD has evolved from 128-bit SSE [16] to 256-bit AVX2 [43] and 512-bit AVX-512 [1]. A single AVX-512 instruction processes eight 64-bit operands in parallel; modern CPU cores provide multiple dedicated SIMD execution units [44], delivering substantially higher data-parallel throughput. The workloads, where SIMD has paid off, linear algebra, computer vision, AI/ML, and modular arithmetic [5, 12, 21, 35, 45, 64, 65], share one property: their sub-tasks are largely independent and map cleanly onto vector lanes1 . Large-number addition, subtraction, and multiplication do not. Their textbook formulations carry strong sequential dependencies between limbs; these dependencies, not SIMD compute throughput, bound the achievable lane parallelism. Addition and Subtraction. The textbook limb-by-limb addition computes 𝑆𝑖 = 𝐴𝑖 + 𝐵𝑖 + 𝐶𝑖 −1 [15, 52], where 𝐴𝑖 and 𝐵𝑖 are the 𝑖-th limbs of the inputs, 𝐶𝑖 −1 is the carry-in from the previous limb, and 𝑆𝑖 is the resulting limb. The carry-out 𝐶𝑖 is generated based on whether the sum exceeds the limb base, thereby creating a carry chain where the carry-out at position 𝑖 becomes the carry-in at position 𝑖+1. Thus, an 𝑚-limb addition reduces to a chain whose latency grows linearly with 𝑚. On the scalar pipeline, the chain hides inside hardware add-with-carry instructions (ADC/SBB [4, 46]) and serves as the basis for the hand-tuned assembly that ships with GMP and OpenSSL. SIMD lanes, on the other hand, are isolated: AVX-512 does not propagate carries across lanes natively [46, 78], and ARM SVE2’s ADCLB/ADCLT [6] only propagate a carry between adjacent even–odd element pairs within a single instruction, so a full 𝑚-limb carry chain still requires a sequence of dependent operations. A direct SIMD port of the carry loop, therefore, has to rebuild the hardware carry chain in software, using compares, masks, and conditional adds. In our measurements (Table 1), the naive SIMD approach has a carry-to-add overhead ratio of 52.1: for every cycle spent on actual addition, over fifty cycles go to sequential carry propagation (full phase-wise breakdown in Table 1). Subtraction shows the same pattern with borrows. Multiplication. Multiplying two 𝑚-limb numbers needs 𝑚 2 word-level products 𝑎𝑖 × 𝑏 𝑗 that have to be accumulated into a 2𝑚-limb result. The two textbook organizations, row-wise schoolbook, and column-wise Comba [15, 18, 52], differ in how they fold these products together, but both run the inner routine through a single accumulator. For larger operands,

Large numbers, a.k.a big numbers, refer to values beyond the standard 64-bit capabilities of common programming languages. These numbers play a fundamental role in scientific computing [7–9], cryptography [27, 31, 38, 49, 61, 66, 70], and various mathematical software tools [60, 71, 81]. Applications in these domains deal with operands ranging from a few hundred bits to thousands of bits [10, 22, 61, 63, 66, 75]. As standard hardware cannot support operands wider than 64 bits [39], software-based approaches such as GNU Multiple Precision Arithmetic Library (GMP) [34], GNU Multiple Precision Floating-Point Reliable library (MPFR) [68], OpenSSL [66], Fast Library for Number Theory (Flint) [28], mpmath [48], Apfloat [76], BigInt [23], SageMath [71] and Mathematica [81] represent large numbers as arrays of fixedsize limbs. A limb is typically a 32-bit or 64-bit unsigned word aligned with the machine width. Two radix styles are common: a saturated representation uses radix 232 or 264

1 Lanes are the independent 128-bit segments that compose SIMD registers (1 for SSE, 2 for AVX2, and 4 for AVX-512), each executing packed data elements in parallel without cross-lane dependencies.

Leveraging SIMD for Accelerating Large-number Arithmetic

Naive SIMD Add Operations

Ren et al. [69] %

Load Add Initial carry detect Sequential carry prop Store

2.6 1.8 1.0 92.8 1.8

Carry/Add ratio

52.1

Operations Load Add Generate carry Add carry Store

Two-level KSA [82] % 5.8 6.7 69.5 13.5 4.5

Operations

Load Add Generate carry Add carry Store

12.4

8.0 8.5 43.8 33.6 6.1 9.1

Naive SIMD Mul Operations Radix conversion Load & permute IFMA compute Store & normalize

Gueron et al. [37] (IFMA) % 8.9 12.7 73.8 4.6 –

Operations Radix conversion Load & permute IFMA compute Store & normalize

% 6.7 9.7 36.7 46.9 –

Table 1. Operation-wise cycles breakdown (%) of prior SIMD arithmetic routines: 512-bit addition and 256-bit multiplication. For addition

routines, the carry-to-add overhead ratio (carry-handling cycles / addition cycles) is also shown. Carry-handling includes: initial carry detection, carry generation, sequential carry propagation, and add carry steps.

libraries fall back to Karatsuba (Algorithm 4 in Appendix A), Toom-Cook, or FFT-based methods [50, 72, 77], but the recursion eventually bottoms out at this same base-case routine. Crucially, Karatsuba’s identity for multiplying two limbs, 𝐴 and 𝐵: 𝐴 × 𝐵 = (𝐴𝐻 × 𝐵𝐻 ) × 𝛽 2 + (𝐴𝐻 +𝐴𝐿 ) × (𝐵𝐻 +𝐵𝐿 ) − 𝐴𝐻 × 𝐵𝐻 − 𝐴𝐿 × 𝐵𝐿 × 𝛽 + 𝐴𝐿 × 𝐵𝐿 , where 𝛽 = 2𝑘 is the limb base, trades one full-size multiplication for three halfsize ones plus additions and subtractions at every level. As recursion deepens, the aggregate addition/subtraction work grows rapidly; faster add/sub therefore improves not only the primitives themselves but also the recursive multiplication, and through it, division, modular exponentiation, and computations like 𝜋. A direct SIMD port of the schoolbook inherits the innerloop dependency: each iteration broadcasts a 𝑏 𝑗 , performs a Fused Multiply-Add (FMA) on the resulting row into the accumulator, and has to finish that fold before the next iteration can start. The routine is not throughput-bound on the multiplier, but latency-bound on the long RAW chain through the accumulator, which leaves the multiple available CPU FMA execution ports underutilized [46](detailed breakdown in Table 1). Appendix A covers the scalar hardware and library context (ALU width, hardware add-with-carry chains, base-case routines in GMP and OpenSSL, and recursive multiplication thresholds) in detail. 2.3

match the carry behavior of the original operands, enabling fast determination of carries across all limbs. The mechanism is clever, but preparing those packed states needs a substantial amount of SIMD-to-scalar conversion. The resulting carry-to-add overhead ratio is 12.4 (Table 1), meaning carry management still dominates runtime by over an order of magnitude. While the authors report gains on isogeny workloads, they do not directly compare their approach against GMP or OpenSSL, while acknowledging the high overheads. Two-level Kogge-Stone. The tool, y-cruncher [83], used for computing the value of mathematical constants like 𝜋 takes a different route: a two-level Kogge-Stone-style addition [53]. While the original KSA is a recursive radix-2 approach with a depth of log 𝑛 for 𝑛-bit operands, they have adopted a twolevel approach with radices of 𝑛𝑘 and 𝑘. In this approach, the operands are divided into 𝑘 groups, each of size 𝑛𝑘 bits, and processed independently at the first level. The second level aggregates results across the 𝑘 groups to compute final carry-bits and adjust sums using carry and max-sum information. The second level sidesteps the per-limb serialization of the naive SIMD path; however, the second-level resolution (carry preparation, mask manipulation, and SIMD-to-scalar transitions) now dominates the runtime. The carry-to-add overhead ratio is 9.1 (Table 1), which is better than naive or Ren et al., but the carry chain’s cost has migrated rather than disappearing. FMA schoolbook with a shared accumulator. For multiplication, prior SIMD work ranges from early SSE2 routines [19] to AVX-512IFMA2 implementations [25, 26, 37, 51]. The state-of-the-art IFMA routine by Gueron et al. [37] adapts the schoolbook organization to AVX-512IFMA: each iteration broadcasts a 𝑏 𝑗 limb, FMAs the resulting row into a single shared accumulator vector, and drains one output limb per iteration through alignr_epi64. The scalar carry chain is gone: all the work now stays in AVX-512 registers. But the accumulator is shared across iterations, creating a long RAW

Prior SIMD-based works

Prior work attempted to break the dependency we discussed above. However, we observe that across all solutions, the sequential chain is removed in one place only to reappear, in a different form, somewhere else: usually at the SIMD-scalar boundary or inside a shared accumulator. Carry-select. Ren et al. [69] target cross-lane carry propagation by simulating carry-select adders [11] to mitigate the long carry dependency chain. Their method breaks down the addition of large integers into smaller, parallel additions of 8-bit operands. Each limb’s addition falls into one of three cases: no-carry, propagate, or generate, determined solely by the smaller operands. The smaller operands are arranged to

2 AVX-512IFMA (Integer Fused Multiply-Add) provides mainly two in-

structions: vpmadd52luq and vpmadd52huq, for performing integer fused multiply-add on 52-bit operands. 3

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

Definition

Algorithm 1: DoT Addition

𝑘

Bit-width of a single limb

𝑛

Total bit length of a large integer operand

𝑚

Total number of limbs needed to represent an operand

𝑤

SIMD vector width (no. of limbs processed in parallel)

Data: Arrays 𝐴, 𝐵 each with 𝑚 limbs of 𝑘 bits Result: Sum array 𝑆 and final carry 𝑐𝑜𝑢𝑡 1 // 𝑛-bit addition: process limbs in chunks of 𝑤 2 Function: DoT-Add-Words(𝐴, 𝐵, 𝑚) 3 begin 4 𝑐𝑜𝑢𝑡 ← 0; 5 for 𝑖 ← 0 to 𝑚 − 1 step 𝑤 do 6 (r, 𝑐𝑜𝑢𝑡 ) ← Add-W-Limbs(𝐴, 𝐵, 𝑖, 𝑐𝑜𝑢𝑡 , 𝑤 ); 7 // Store Results 8 simd_store(𝑤 limbs of r into 𝑆 [𝑖 ] to 𝑆 [𝑖 + 𝑤 − 1] ); 9 end 10 return (𝑆, 𝑐𝑜𝑢𝑡 ); 11 end 12 // Function for adding 1 ≤ 𝑤 ≤ 8 limbs simultaneously 13 Function: Add-W-Limbs(𝐴, 𝐵, 𝑖, 𝑐 𝑖𝑛 , 𝑤 ) 14 begin 15 max ← vector of 𝑤 limbs, each with value 2𝑘 − 1; 16 // Phase 1: Parallel Addition 17 a ← simd_load(𝑤 limbs from 𝐴[𝑖 ] to 𝐴[𝑖 + 𝑤 − 1] ); 18 b ← simd_load(𝑤 limbs from 𝐵 [𝑖 ] to 𝐵 [𝑖 + 𝑤 − 1] ); 19 r ← simd_add(a, b) ; 20 // Phase 2: Carry Generation 21 c ← simd_cmp_lt(r, a) ; // Generate mask 22 𝑐𝑜𝑢𝑡 ← c ≫ (𝑤 − 1); 23 c ← (c ≪ 1) | 𝑐𝑖𝑛 ; 24 // Phase 3: Carry Propagation 25 r′ ← simd_add_carry(r, c); 26 c ← simd_cmp_lt(r′ , r) ; // Generate mask 27 if c ≠ 0 then 28 // Phase 4: Carry Adjustment 29 c ← c ≪ 1; 30 m ← simd_cmp_eq(r′ , max) ; // Generate mask 31 c ← c + m; 32 𝑐𝑜𝑢𝑡 ← 𝑐𝑜𝑢𝑡 | (c ≫ 𝑤 ); 33 m ← c ⊕ m; 34 r′ ← simd_add_carry(r′ , m); 35 end 36 return (r′ , 𝑐𝑜𝑢𝑡 ); 37 end

Symbol

a, b, r, ...

SIMD vector registers

𝐴𝑖 , 𝐵𝑖

𝑘-bit values

Table 2. Notation used throughout the DoT algorithm descriptions. Boldface lowercase letters (a, b, r, . . . ) denote SIMD vector registers; uppercase subscripted letters (𝐴𝑖 , 𝐵𝑖 ) denote individual 𝑘-bit limb values.

dependency chain through madd52 → alignr → madd52. Although newer x86-64 CPUs expose two IFMA execution ports [44], only one can issue per cycle along this chain, so the routine is latency-bound rather than throughput-bound, achieving an IPC of only 4.2 despite the availability of two IFMA ports. A full phase-wise breakdown is in Table 1. We could not directly evaluate the other IFMA-based works [24– 26] due to the lack of source code and insufficient detail in the papers to reimplement them faithfully. Key Takeaway. Across all methods in Table 1, the SIMD compute phase is never what limits performance. The surrounding work needed to feed the intermediate computations around a dependency imposed by the sequential algorithm dominates the computation. None of the prior work goes back and restructures the operation around better independent, data-parallel work in the first place. For addition, they have tried to eliminate the sequential carry dependency, either by paying a heavy SIMD-to-scalar transition for preparing the carry states or by having a complex two-level reduction, diminishing the SIMD benefits. On the multiplication side, the dependency through the shared accumulator creates a long RAW chain that starves the available CPU execution ports. Two natural questions follow: (1) Are carry cascades common enough in practice to justify paying the high overhead of complex carry-handling costs? (2) Can multiplication’s partial products be reorganized to eliminate the sharedaccumulator dependency?

while guaranteeing correctness under all inputs and microarchitecture portability. 3.1

DoT Addition and Subtraction

As discussed earlier, addition requires carry propagation, which introduces necessary dependencies and limits the SIMD utilization. Carefully analyzing the random nature of large number inputs and the limb representation, we observe that propagating carry-chains is rare in practice; a fact not accounted for by prior work. Using a limb-based approach, addition can be limited to a single neighboring carry-propagation. In other words, corresponding limbs of both operands can be added in parallel, followed by the generated carry-bits to be propagated to preceding intermediate sums in parallel, too. In the majority of

Design

We describe the design of DigitsOnTurbo (DoT), an SIMDbased algorithm for large-number addition (and subtraction3 ), and multiplication. In both cases, the design starts by restructuring the computation around independent, dataparallel operations, rather than vectorizing an existing sequential implementation, and builds SIMD phases around that structure. Our design focuses on maximizing parallelism 3Without loss of generality, we focus on addition in this paper; differences

for subtraction are noted where they arise.

Leveraging SIMD for Accelerating Large-number Arithmetic

Phase 3. The aligned carries are added to the intermediate sums in parallel: 𝑅𝑖′ = 𝑅𝑖 + 𝑐𝑖 (line 25). The comparison on line 26 checks whether any of the carry additions overflow. If not (the common case), the result is stored and the function returns. Otherwise, Phase 4 handles the rare cascade. Phase 4. The secondary carry-bits are shifted left (line 29) and combined with a mask of limbs equal to 2𝑘 −1 (lines 30– 33), using the carry-adjustment trick from the Kogge-Stone adder [53]. The adjusted carry is added back into the sums (line 34) to give the final result. Because Phase 4 operates on values that have already been carry-propagated in Phase 3; the adjustment is bounded and correct. When 𝑚 is not a multiple of 𝑤, the final call to AddW-Limbs uses masked SIMD load/store variants to process only the remaining limbs; smaller operands can also use a narrower 𝑤 (e.g., 𝑤=4 for AVX2 or 𝑤=2 for SSE) by selecting the corresponding intrinsic width.

Figure 1. Illustration of DoT addition for a 4-limb example. Phase 1 (P1) and Phase 3 (P3) perform SIMD ADD in parallel; Phase 2 (P2) generates and shifts carry-bits on scalar/mask registers; Phase 4 (P4) handles the rare carry-cascade case via the slow path.

cases, propagating carry-bit to preceding intermediate sums may not generate an additional carry-bit. A new carry-bit is generated only when the earlier carry of 1 is added to an intermediate sum that equals the maximum value for the base (e.g., 264 − 1 for 64-bit limbs). Such propagation cascades only when all intermediate limbs compute to the maximum base-value. However, the probability of each of the sums computing to the maximum possible value is quite less ( 21𝑘 where 𝑘 is the number of bits in a limb), which allows us to avoid these cascading propagations in majority of the cases and defer them to a slow path that only executes when such a cascade occurs. Appendix B provides a formal analysis of the same, showing that the probability of a carry cascade is negligible for random inputs. Considering the above observation, we propose a fourphase DoT Addition algorithm. Phase 1 (P1) performs limbwise addition of 𝐴𝑖 and 𝐵𝑖 in parallel, producing intermediate sums without any carry management. Phase 2 (P2) detects the carries from P1 and shifts them one position to align with the limbs they must propagate into, with the top-limb carry extracted as 𝑐𝑜𝑢𝑡 . Phase 3 (P3) adds these aligned carries (and any incoming 𝑐𝑖𝑛 ) to the intermediate sums in a single parallel step. Phase 4 (P4) handles the cascading case where P3 itself generates a new carry by adjusting the sums with a carry mask. Figure 1 illustrates the flow for a 4-limb example. Algorithm 1 shows the pseudocode. We walk through the algorithm with 𝑛=512-bit operands and 𝑘=64-bit limbs, so 𝑚=8 limbs fit in a single AVX-512 call with 𝑤=8. DoT-AddWords calls Add-W-Limbs for 𝐴0 –𝐴7 , 𝐵 0 –𝐵 7 , with 𝑐𝑖𝑛 = 0. Phase 1. All eight limbs of 𝐴 and 𝐵 are loaded into SIMD registers and added in parallel (lines 17–19), computing 𝑅𝑖 = 𝐴𝑖 + 𝐵𝑖 for 𝑖 = 0..7 in one instruction. Deferring carry detection ensures no inter-lane dependency. Phase 2. A carry at position 𝑖 is detected by 𝑐𝑖 = (𝑅𝑖 < 𝐴𝑖 ) via a SIMD compare across all eight lanes (line 21). The carry out of the top limb (𝑐 7 ) is saved as 𝑐𝑜𝑢𝑡 (line 22); the remaining bits are shifted left by one and the incoming 𝑐𝑖𝑛 is inserted at bit 0 (line 23), giving c = 𝑐 6 . . . 𝑐 0 𝑐𝑖𝑛 and aligning each carry-bit with the limb it must propagate into.

How DoT Reduces Carry Overhead. Prior SIMD-based solutions do not handle common and rare cases separately; consequently, they pay heavy overhead for carry management even when no cascade occurs. On the contrary, DoT’s Phase 2 performs three simple tasks: detect the carry-bits using a SIMD compare, shift the carry by one position to align with the limbs they must propagate into, and extract the top-limb carry (𝑐𝑜𝑢𝑡 ). For common cases, only the first three phases are executed. Phase 4 is executed rarely when a carry cascade arises. By isolating the carry propagation to a separate phase that only executes when necessary, DoT significantly improves common case performance. Proof of correctness. We show that Dot-Add-Words is correct and produces the same result as normal addition with carry across limbs. Theorem 3.1 states the correctness theorem; the proof is provided in the Appendix. Theorem 3.1 (Correctness of DoT-Addition). Suppose two Í 𝑖 𝑚−1 + · · · + large integers 𝐴 and 𝐵, 𝐴 = 𝑚−1 𝑖=0 𝐴𝑖 𝑋 = 𝐴𝑚−1𝑋 Í 𝑚−1 1 0 𝑗 𝑛−1 𝐴1𝑋 + 𝐴0𝑋 , 𝐵 = 𝑗=0 𝐵 𝑗 𝑋 = 𝐵𝑚−1𝑋 + · · · + 𝐵 1𝑋 1 + 𝑘 0 𝐵 0𝑋 , where 𝑋 = 2 (where 𝑘 is the bit size of a limb), the number of limbs for the two integers is 𝑚, and 0 ≤ 𝐴𝑖 , 𝐵 𝑗 < 𝑋 . Í 𝑖 If DoT-Add-Words(𝐴, 𝐵, 𝑚) = (𝑆, 𝑐 out ), then 𝑆 = 𝑚−1 𝑖=0 (𝑟 𝑖 𝑋 ) such that 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋, and 𝑐 out = 𝑐𝑚 . Subtraction. DoT subtraction mirrors addition, replacing carries with borrows. Limbs are subtracted in parallel in Phase 1; borrow bits are generated and aligned in Phase 2; borrows are propagated in Phase 3; and Phase 4 handles the case where a limb that reached zero in Phase 1 receives a borrow in Phase 3, triggering further borrow generation. The probability of Phase 4 is similar to the addition case. 3.2

DoT Multiplication

Large-number multiplication presents a different challenge of computing and accumulating partial products efficiently 5

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

without serialization. To overcome the challenge, we adopt the vertical and crosswise multiplication technique [57, 62, 80], which performs multiplication column-by-column, i.e., for each output position 𝑐, all cross-products 𝐴𝑖 × 𝐵 𝑗 , such that 𝑖 + 𝑗 = 𝑐, are computed and summed together, i.e., for 𝑐 = 0, we have the cross-product 𝐴0𝑋 0 .𝐵 0𝑋 0 , for 𝑐 = 1, we have the cross-product (𝐴1𝑋 1 .𝐵 0𝑋 0 + 𝐴0𝑋 0 .𝐵 1𝑋 1 ) and so on.

Algorithm 2: DoT Multiplication Data: Arrays 𝐴, 𝐵 each with 𝑚 limbs of 𝑘 bits Result: Product array 𝑃 = 𝐴 × 𝐵 Function: DoT-Mul-Words(𝐴, 𝐵, 𝑚, 𝑘 ) begin 3 𝑀𝑎 [0 . . . 𝑚 2 − 1], 𝑀𝑏 [0 . . . 𝑚 2 − 1], 𝑃 _𝑙𝑜 [0 . . . 𝑚 2 − 1], 𝑃 _ℎ𝑖 [0 . . . 𝑚 2 − 1], 𝑐𝑜𝑙_𝑙𝑜 [0 . . . 2𝑚 − 1] ← 0; 4 𝑖𝑑𝑥_𝑎 ← 0, 𝑖𝑑𝑥_𝑏 ← 0; 5 // Phase 1: Gather limbs into columns 6 for 𝑐 ← 0 to 2𝑚 − 2 do 7 for all pairs (𝑖, 𝑗 ) such that 𝑖 + 𝑗 = 𝑐 do 8 𝑀𝑎 [𝑖𝑑𝑥_𝑎 + +] ← 𝐴[𝑖 ]; 9 𝑀𝑏 [𝑖𝑑𝑥_𝑏 + +] ← 𝐵 [ 𝑗 ]; 10 end 11 end 12 // Phase 2: Compute partial products with SIMD 13 for 𝑖 ← 0 to 𝑚 2 − 1 step 𝑤 do 14 a ← simd_load(𝑤 limbs from 𝑀𝑎 [𝑖 ] to 𝑀𝑎 [𝑖 + 𝑤 − 1] ); 15 b ← simd_load(𝑤 limbs from 𝑀𝑏 [𝑖 ] to 𝑀𝑏 [𝑖 + 𝑤 − 1] ); 16 p_lo ← simd_mul_lo(a, b); 17 p_hi ← simd_mul_hi(a, b); 18 simd_store(p_lo into 𝑃 _𝑙𝑜 [𝑖 ] to 𝑃 _𝑙𝑜 [𝑖 + 𝑤 − 1] ); 19 simd_store(p_hi into 𝑃 _ℎ𝑖 [𝑖 ] to 𝑃 _ℎ𝑖 [𝑖 + 𝑤 − 1] ); 20 end 21 // Phase 3: Align hi parts to neighboring columns 22 𝑖𝑑𝑥 ← 0; 23 for 𝑐 ← 0 to 2𝑚 − 2 do 24 𝑛𝑢𝑚_𝑝𝑎𝑖𝑟𝑠 ← min(𝑐 + 1, 𝑚, 2𝑚 − 1 − 𝑐 ); 25 for 𝑝 ← 0 to 𝑛𝑢𝑚_𝑝𝑎𝑖𝑟𝑠 − 1 do 26 𝑡 ← 𝑖𝑑𝑥 + 𝑝; 27 𝑐𝑜𝑙_𝑙𝑜 [𝑐 ] ← 𝑐𝑜𝑙_𝑙𝑜 [𝑐 ] + 𝑃 _𝑙𝑜 [𝑡 ]; 28 𝑐𝑜𝑙_𝑙𝑜 [𝑐 + 1] ← 𝑐𝑜𝑙_𝑙𝑜 [𝑐 + 1] + 𝑃 _ℎ𝑖 [𝑡 ]; 29 end 30 𝑖𝑑𝑥 ← 𝑖𝑑𝑥 + 𝑛𝑢𝑚_𝑝𝑎𝑖𝑟𝑠; 31 end 32 // Phase 4: Reduce partial products in each column 33 for 𝑐 ← 0 to 2𝑚 − 1 do 34 𝑃 [𝑐 ] ← 𝑐𝑜𝑙_𝑙𝑜 [𝑐 ]; 35 end 36 // Phase 5: Carry-over adjustment 37 carry ← 0; 38 for 𝑐 ← 0 to 2𝑚 − 1 do 39 𝑃 [𝑐 ] ← 𝑃 [𝑐 ] + carry; 40 carry ← 𝑃 [𝑐 ] ≫ 𝑘; 41 𝑃 [𝑐 ] ← 𝑃 [𝑐 ] mod 2𝑘 ; 42 end 43 if carry > 0 then 44 Prepend carry as an extra limb to the result; 45 end 46 return 𝑃 ; 47 end 1 2

Figure 2. “Vertical and Crosswise” partial product organization

for 2×2, 3×3, and 5×5 limb multiplication. Each line represents one cross-product 𝐴𝑖 × 𝐵 𝑗 ; lines of the same color belong to output column 𝑐 = 𝑖+𝑗 and are summed together. A 2𝑚−1-column structure exposes all 𝑚 2 partial products as independent computations.

Crucially, all cross-products are independent of one another, so they can all be computed before any summation or carry adjustment. Figure 2 illustrates this for 2×2, 3×3, and 5×5 limb multiplications. A final carry-adjustment pass propagates overflows from each column into the next to produce the final product. Since all partial products are independent, 𝑤 limb pairs can be gathered into SIMD registers per output column and computed in a single instruction. A horizontal reduction sums the partial products per column, followed by the carry-adjustment pass. The vertical and crosswise multiplication is defined as: VnC(𝐴, 𝐵, 𝑛, 𝑘) = 𝐴0 ·𝐵 0 + (𝐴1 2𝑘 ·𝐵 0 + 𝐴0 ·𝐵 1 2𝑘 ) + · · · + (𝐴𝑛−1 2𝑘 (𝑛−1) ·𝐵𝑛−2 2𝑘 (𝑛−2) + 𝐴𝑛−2 2𝑘 (𝑛−2) ·𝐵𝑛−1 2𝑘 (𝑛−1) ) + (𝐴𝑛−1 2𝑘 (𝑛−1) ·𝐵𝑛−1 2𝑘 (𝑛−1) ). DoT Multiplication algorithm comprises five phases. The pseudocode is provided in Algorithm 2. Phase 1 gathers the limb pairs for each output column 𝑐 = 𝑖 + 𝑗 into flat arrays 𝑀𝑎 and 𝑀𝑏 . Phase 2 computes every partial product in parallel using SIMD vector multiplies with a zero accumulator, eliminating the serialized FMA chain of prior IFMA routines [37] and fully utilizing both IFMA execution ports. Phase 3 propagates each high half to the next column’s accumulator via a single flat-index traversal. Phase 4 performs a horizontal reduction per column, and Phase 5 executes the carry-adjustment pass, which performs a single sequential pass over the 2𝑚 − 1 column sums. We walk through Algorithm 2 with 𝑚=5 and 𝑘=52: two 260-bit operands 𝐴=𝐴4 . . . 𝐴0 and 𝐵=𝐵 4 . . . 𝐵 0 produce a 520bit product 𝑃 across 2𝑚−1=9 columns with 1, 2, 3, 4, 5, 4, 3, 2, 1 pairs per column respectively (the triangular pattern in Figure 2). DoT-Mul-Words is invoked on 𝐴0 –𝐴4 , 𝐵 0 –𝐵 4 and returns a 10-limb product.

Phase 1. The gather loop (lines 6–9) enumerates each output column 𝑐 and every (𝑖, 𝑗) with 𝑖+𝑗=𝑐, packing operands into flat arrays 𝑀𝑎 and 𝑀𝑏 in column-wise groups. For 𝑐=4, it emits five pairs: (𝐴0, 𝐵 4 ), (𝐴1, 𝐵 3 ), (𝐴2, 𝐵 2 ), (𝐴3, 𝐵 1 ), (𝐴4, 𝐵 0 ); after the loop, 𝑀𝑎 and 𝑀𝑏 each hold 25 entries with every 6

Leveraging SIMD for Accelerating Large-number Arithmetic

DoT Addition (𝑚 = 8, 64-bit limbs) Operations

DoT Mul (𝑚 = 4, 64-bit limbs)

Random %

Pathological %

Load (Phase 1) Add (Phase 1) Generate carry (Phase 2) Add carry (Phase 3) Store & check (Phase 3) Overflow handling (Phase 4)

18.7 13.8 25.8 20.7 21.0 —

7.3 5.4 10.1 8.1 8.2 60.9

Carry/Add ratio

4.9

16.2

DoT Mul (𝑚 = 5, 52-bit limbs)

Operations

Gather & radix conversion (Phase 1) IFMA compute (Phase 2) Align hi parts (Phase 3) Reduce columns (Phase 4) Carry-over & store (Phase 5)

22.6 19.4 23.5 24.4 10.1

Operations Gather (Phase 1) IFMA compute (Phase 2) Align hi parts (Phase 3) Reduce columns (Phase 4) Carry-over & store (Phase 5)

—

% 8.0 22.9 27.3 31.5 10.3 —

Table 3. Phase-wise timing breakdown (%) of DoT for 512-bit addition (𝑚=8, 64-bit limbs) and 256-bit multiplication (𝑚=4 and 𝑚=5 limbs). For addition, the random column shows percentages for random inputs (Phase 4 never fires); the pathological column shows percentages for pathological inputs (Phase 4 fires).

Í 𝑗 𝑚−1 + 𝐴𝑚−1𝑋 𝑚−1 +· · ·+𝐴1𝑋 1 +𝐴0𝑋 0 , 𝐵 = 𝑚−1 𝑗=0 𝐵 𝑗 𝑋 = 𝐵𝑚−1𝑋 · · · + 𝐵 1𝑋 1 + 𝐵 0𝑋 0 , where 𝑋 = 2𝑘 (where 𝑘 is the bit size of a limb), the number of limbs for the two integers is 𝑚, and 0 ≤ 𝐴𝑖 , 𝐵 𝑗 < 𝑋 . Then, 𝐴 × 𝐵 = VnC(𝐴, 𝐵).

cross-product exposed as an independent operand pair with no dependencies between them. Intuitively, by restructuring the computation around the output columns, all pairs contributing to a column can be computed in parallel in Phase 2 directly without shuffling or permutating. Phase 2. The 25 entries are consumed in chunks of 𝑤=8 via SIMD multiplies (lines 14–19): three full iterations cover 24 pairs plus one trailing pair handled as the remainder step. Each vector issues simd_mul_lo/simd_mul_hi against a zero accumulator, producing eight independent low halves and eight high halves per call. Because no iteration reads another’s accumulator, both IFMA ports stay busy in parallel with no RAW chain through the sequence, unlike Gueron et al.’s shared-accumulator chain [37]. Phase 3. The alignment pass (lines 23–28) folds each product into its column: 𝑃_𝑙𝑜 [𝑡] is added into 𝑐𝑜𝑙_𝑙𝑜 [𝑐] while 𝑃_ℎ𝑖 [𝑡] is promoted into 𝑐𝑜𝑙_𝑙𝑜 [𝑐+1]. The pair (2, 3), for instance, deposits its low half at column 5 and its high half at column 6. This phase is a simple flat loop with no dependencies between iterations, and thus can be parallelized with SIMD or left as a scalar pass. The result is that each column 𝑐 holds the full sum of its cross-products, albeit with potential overflow above 2𝑘 that must be handled in Phase 5. Phase 4. Each column’s accumulated value is moved into the output array 𝑃 [𝑐] (line 34). Phase 5. A single scalar pass over the 2𝑚 − 1=9 column sums (lines 38–41) propagates each column’s overflow above 252 into the next, yielding the final 9-limb product and an optional overflow limb (line 44). Thus, only Phase 5 is sequential; Phases 1–4 are fully data-parallel and costlier serial work is limited to this short tail pass.

3.3

Implementation

We implement DoT in C using x86-64 SIMD intrinsics [47], trading peak microarchitectural tuning for portability across compilers and microarchitectures. Addition and subtraction support saturated radix (𝑘 = 64, base 264 ), using masked load/store intrinsics to handle non-multiple-of-𝑤 limb counts, and can be adopted for SSE (𝑤 = 2), AVX2 (𝑤 = 4), and AVX512 (𝑤 = 8) SIMD widths by selecting the corresponding intrinsics. Multiplication uses two fixed-size implementations, 5×5 and 4×4, both realizing Algorithm 2’s five phases. The 5×5 routine operates on unsaturated radix (𝑘 = 52, base 252 ), since AVX-512 IFMA instructions operate on 52-bit limbs. The 4×4 routine is designed for compatibility with the saturated radix (𝑘 = 64, base 264 ) used by GMP and OpenSSL, and pays the extra cost of radix conversion packing at entry (𝑘 = 64 → 𝑘 = 52) and unpacking (𝑘 = 52 → 𝑘 = 64) at exit. We integrate DoT as a static library into the C pathways of GMP-6.3.0 (DoTMP) and OpenSSL-3.5.2 (DoTSSL), replacing each library’s add/sub primitives and base-case multiplication; higher-level recursive multiplications automatically invoke DoT as the base case. Unaligned SIMD load/store intrinsics throughout preserve compatibility with both libraries’ non-aligned memory layouts. We also test DoT on the integrated libraries’ test suites, which cover a wide range of operand sizes and values, and a hefty set of higher-level applications (e.g., RSA, DSA) that rely on the primitives. In total, DoT consists of 1013 lines of C code, with 38 and 43 lines of integration changes in GMP and OpenSSL, respectively4 .

Proof of correctness. The correctness of vertical and crosswise multiplication is shown as Theorem 3.2. Algorithm 2 is an adaptation of the vertical and crosswise multiplication approach designed to work with the SIMD architecture and the available registers. The proof is shown in the Appendix. Theorem 3.2 (Correctness of vertical and crosswise multipliÍ 𝑖 cation). Suppose two large integers 𝐴 and 𝐵, 𝐴 = 𝑚−1 𝑖=0 𝐴𝑖 𝑋 =

4 Our implementation, library variants, and experimental setup are available

at: https://anonymous.4open.science/r/DigitsOnTurbo/.

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

release their implementation, we derived a fair and AVX512-intrinsics-based implementation of their ProposedAdd algorithm from the pseudocode provided in their paper.

Evaluation

We evaluate DoT along four axes: (i) reduction in carrypreparation overhead for add/sub, (ii) scaling across SIMD widths (𝑤=2, 4, 8) for add/sub, (iii) performance of base-case multiplication, (iv) impact on higher-level operations. We run the experiments on an Intel Xeon Gold 6548Y+ CPU (Emerald Rapids microarchitecture) with 64 cores (two sockets, 32 cores per socket) @ 2.50 GHz, 256 GB of DRAM, and Ubuntu 22.04.2 LTS. We also run experiments on an Intel Xeon Max 9462 (Sapphire Rapids, 64 cores @ 2.70 GHz, same OS and DRAM). We report results from the 6548Y+ in the main paper and include results from the 9462 in the Appendix. We use clang version 14.0.0-1ubuntu1.1 with −𝑂2 and −𝑚𝑎𝑟𝑐ℎ = 𝑛𝑎𝑡𝑖𝑣𝑒 optimization flags to compile DoT. To evaluate the accuracy and performance of DoT, we employ two sets of test cases: random and pathological. Random cases mimic typical usage; pathological cases target edge scenarios (full carry/borrow propagation, maxed-out/zero limbs, frequent carries, frequent borrows, and mixed cases). We use the Mersenne Twister [59] seeded with random integers to generate the test cases. Every test run includes 100,000 random and 1,000 pathological test cases for twelve operand sizes from 512 to 32768 bits. Micro-benchmark results are reported as the mean over twenty runs per experimental setup (a unique combination of operation and operand size), using the 95% confidence interval. We compute timing and throughput based on GMPbench’s methodology [33] and CPU ticks via the RDTSC instruction [67]. Speedups, wherever reported, are measured as the ratio of execution times, which are also verified using CPU ticks and throughputs. Instruction counts, wherever reported, were measured using the perf_event_open system call [58]. Unless otherwise noted, all compared methods are evaluated in the same benchmarking harness with identical operand sets, iteration counts, and measurement methodology (RDTSC and perf_event_open); each implementation is validated against precomputed expected outputs.

4.1

Implementation dot_mul_5 × 5 dot_mul_4 × 4 OpenSSL BN_mul Gueron & Krasnov [37] GMP mpz_mul

Instructions

Avg. Cycles

IPC

221 265 655 345 870

28.6 35.2 68.7 81.5 88.4

7.7 7.5 9.5 4.2 9.8

Table 4. Instruction counts, average cycle counts, and IPC for 256-bit multiplication. Cycle counts were obtained via RDTSC

Randomly Generated Test Cases. Compared to the TwoLevel KSA, DoT addition achieves a speedup of 1× to 1.9×. For subtraction, DoT achieves a speedup of 0.9× to 1.9×. DoT, on average, reduces instruction count by 15% for addition and 16% for subtraction. Compared to Ren et al., DoT achieves a speedup of 1.4× to 2.2× and 1.1× to 1.8× for addition and subtraction, respectively. DoT reduces the instruction count by an average of 28.8% for addition and 17.3% for subtraction, respectively. For addition and subtraction, DoT shows modest speedups at small operand sizes due to fixed function-call overhead. As operands grow, the overhead amortizes, and the lower per-limb computation cost dominates. Pathological Test Cases. Due to space constraints, the plot comparing the performance for pathological test cases is shown in the Appendix C. DoT achieves a speedup 0.7× to 2× against Two-level KSA and 0.8× to 1.9× against Ren et al. 4.2

DoT’s Performance over SIMD Widths

We now compare the performance of three DoT SIMD variants (𝑤 = 2, 4, and 8) with the scalar non-SIMD variant for addition and subtraction. Figure 3(b) shows these results. For comparing the baseline scalar, we have used _addcarryx _u64 and _subborrow_u64 intrinsics. With SSE (𝑤 = 2), we observe a geometric mean speedup of 0.7× for both addition and subtraction; at lower SIMD width, the overhead outweighs the parallelism benefits. With AVX2 (𝑤 = 4), the geometric mean speedups increase to 1.2× for addition and subtraction, respectively. The AVX512 (𝑤 = 8) variant achieves the highest speedups, with geometric mean speedups of 1.8× for addition and subtraction, respectively. As stated earlier, for the lower range of the operands (512– 4096 bits), the speedups are more modest, while for larger operands (6144–32768 bits), the performance gains are more significant and consistent. SSE and AVX2 increase instruction count by 125% and 25% for addition and 129% and 29% for subtraction over scalar, respectively. AVX512 reduces instruction count by 25% for addition and 22% for subtraction, reaching 33% and 30% reductions for the largest operands. These results reveal that SIMD parallelism and not instruction compactness drives

Reduction of Carry Preparation Overhead

Table 1 and Table 3 give the full phase-wise breakdown for prior approaches and DoT, respectively. For random inputs (Phase 4 never triggered across our entire test set), DoT achieves a carry-to-add ratio of ∼ 4.9, roughly half that of two-level KSA (∼ 9.1) and a third of Ren et al. (∼ 12.4). DoT also runs 1.85× faster overall, so the absolute carry cycles per limb are even lower. For pathological inputs where Phase 4 fires, the ratio rises to 16.2, still well below naive SIMD (52.1). Figure 3(a) shows the performance comparison of DoT (AVX512), two-level KSA adapted from y-cruncher’s approach, and Ren et al.’s method [69] for addition and subtraction on both randomly generated test cases. Since Ren et al. did not 8

Leveraging SIMD for Accelerating Large-number Arithmetic

Method DoT KSA Ren et al.

8 76

Mul (DoTSSL/OpenSSL)

Bit Size

(d) Mul in integrated libs

6 24

8 28 12

96 40

2 51

8 76

6 57

57 24

2 51

Speedup

1.2x

38 16

Mul (DoTMP/GMP)

44 61

72 30

1.4x

1.0x

0x 15

8 76

Bit Size

(b) Add SIMD variants speedup

Operation Add Sub

Bit Size

1.0x

(a) Add/sub vs prior works

1.5x

0.5x

Pair DoTMP/GMP DoTSSL/OpenSSL

2.0x

92 81

96 40

Operation Add Sub

SSE (w=2) AVX2 (w=4) AVX512 (w=8)

Speedup

100

Time (ns) (log scale)

2.5x

Figure 3. Micro-benchmark evaluation of DoT across four axes. (a) Execution time (log scale) of DoT (AVX512), two-level KSA, and Ren et al.

for add/sub across 512–32768-bit random operands. (b) Execution time speedup of DoT SIMD variants (𝑤=2, 4, 8) over scalar add-with-carry. (c) Execution time speedup of DoTMP over GMP and DoTSSL over OpenSSL for add/sub. (d) Execution time speedup of DoTMP over GMP and DoTSSL over OpenSSL for multiplication.

the speedup. The benefit of AVX-512 DoT is not that it issues fewer instructions, but that each SIMD instruction operates on 𝑤 = 8 limbs at once. SSE and AVX2 do not have enough lanes to amortize the carry-management overhead; AVX512 (𝑤 = 8) crosses the threshold where lane parallelism dominates the overhead.

4.3

Despite being faster, DoT spends significantly less time in the IFMA compute phase (19.4%) compared to Gueron and Krasnov’s method (36.7%) (Table 3): DoT’s reorganization of the multiplication routine allows it to better utilize the CPU’s execution resources and reduce dependency chains, leading to higher IPC and overall performance gains. While Gueron and Krasnov target their multiplication routine for 1024-bit and above operand sizes, we included their method for comparison since it is the only publicly available AVX-512-based multiplication implementation; other implementations [24– 26] were not publicly available and their papers did not provide sufficient details to reimplement their algorithms for a fair comparison. Notably, all of the prior works on SIMD multiplication keep GMP as the baseline for comparison, and also observe speedups only beyond 1024-bit operands.

256-bit Multiplication with DoT

Table 4 compares our 256-bit base case multiplication routines (4 × 4 and 5 × 5) against GMP, OpenSSL, and Gueron and Krasnov’s AVX-512IFMA implementation [37]. DoT’s 4 × 4 routine issues 23% fewer instructions than Gueron and Krasnov’s achieving a 2.31× speedup, driven primarily by lower cycle count and higher IPC: 7.5 for DoT compared to 4.2 of Gueron and Krasnov. Against OpenSSL and GMP, the same routine delivers speedups of 1.95× and 2.51×, respectively. Since DoT’s 4 × 4 routine converts the 64-bit limbs to a 52-bit representation to better fit the IFMA instruction’s 52-bit operand limit, it incurs some overhead for the radix conversion and alignment of the hi parts, which is reflected in the timing breakdown in Table 3. The primary reason for implementing the 4 × 4 routine is to have a direct APIcompatibility for the 256-bit multiplication base cases in GMP and OpenSSL, which use 64-bit limb representation.

4.4

DoT-integrated GMP and OpenSSL

Figure 3(c) shows the timing speedup of DoTMP over GMP and DoTSSL over OpenSSL for addition and subtraction across 512–32768-bit operands. For addition, DoTMP and DoTSSL achieve geometric mean speedups of 3.81× and 2.95×, respectively, with larger gains at larger operand sizes. For subtraction, the geometric mean speedups are 3.73× and 4.08×. 9

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

Instruction counts drop substantially: DoTMP sees reductions of 78.6% and 77.4% for addition and subtraction, while DoTSSL sees 65.9% and 71.5%. For multiplication (Figure 3(d)), DoTMP achieves a geometric mean speedup of 1.41× over GMP across 512–32768-bit operands, reducing instruction count by 47.3% on average (up to 53.7%). DoT only replaces the 256-bit base case ( 2.5× faster for GMP), accounting for 50% of total cycles at large operand sizes (the rest spent in higher-level recursion), bounding the speedup to the 1.34–1.49× range. DoTSSL achieves 1.20× geometric mean speedup over OpenSSL (ranging from 1.10× to 1.32×), with instruction count reductions averaging 41.6%. As we integrate DoT multiplication as the base case for 256-bit, multiplication speedups are larger when operand bit sizes are powers of two and take the path to the optimized 256-bit multiplication in GMP’s Toom-Cook/FFT variants and OpenSSL’s Karatsuba implementation. Beyond 4096 bits, GMP switches to unequal partitioning for Toom-Cook, which results in fewer 256-bit multiplications and therefore smaller speedups from DoT multiplication. OpenSSL’s Karatsuba implementation, on the other hand, continues to use the same base case for all operand sizes, resulting in consistent speedups for DoTSSL multiplication over OpenSSL. 4.5

dot_add_words, dot_sub_words, and dot_mul_4x4 (refer to Figure 6). In GMPbench, DoT’s routines account for a significant portion of the cycles in the multiply, divide, and 𝜋 workloads for larger operand sizes, which explains the larger speedups observed in these workloads. In OpenSSL speed, DoT’s routines contribute a smaller but still meaningful share of cycles across RSA, FFDH, and DSA benchmarks, consistent with the more modest speedups observed there.

Discussion

Futuristic Hardware. Existing SIMD support comes with a few limitations. DoT is designed around these limitations to ensure the parallelism offered by SIMD is fully utilized. (i) The scalar add-with-carry (ADC/SBB) instruction eliminates all the explicit carry-handling steps. However, in SIMD, none of the ISAs have native carry propagation across lanes, so all the approaches, including DoT, manage carries explicitly. ARM SVE2 in particular exposes native carry-generation instructions (ADCLB/ADCLT [6]) that map directly to DoT Add Phase 2, suggesting an SVE2 port would make Phase 2 cheaper than the AVX-512 version. (ii) Currently, the SIMD multiplication instruction sets (e.g., AVX-512F), which work with 64-bit limbs, only provide the lower 64-bit product of the 128-bit result. On the other hand, AVX-512IFMA provides a fused multiply-add instruction that operates on lower 52bit limbs, and requires two separate instructions to get the full 104-bit product. As the hardware evolves and new features are added to SIMD, such as add-with-carry chaining or mulx-like widemultiply instruction that gives the full product in one step, approaches like DoT that restructure algorithms around independent, data-parallel operations would further reduce the overhead and make SIMD usage more attractive for largenumber libraries. As DoT uses C with AVX-512 intrinsics and no microarchitecture-specific assembly, future compiler improvements would benefit DoT without code changes. AVX10 [42], a successor to AVX-512, preserves the intrinsics of AVX-512, DoT will run unchanged on upcoming Intel Pand E-cores.

DoT’s Impact on Higher-Level Applications

GMPbench. Figure 4 shows DoTMP’s improvement across the GMPbench suite. DoTMP improves the overall GMPbench score by 7.8%. The multiply aggregate improves by 15.3% (up to 48.7% for 2, 097, 152 × 2, 097, 152-bit operands). Divide aggregate improves by 8.4% (up to 36.4%) even though DoT replaces no division routine, because GMP’s NewtonRaphson divider calls multiplication and add/sub in tight loops. Computation of 𝜋 gains 13.3% (up to 19.3% for 1M digits) because Chudnovsky leans almost entirely on largeinteger multiply and divide. Even GCD gains by 3.1% (up to 13.7%), since GMP’s Lehmer-Euclid hybrid bottoms out in large addition and subtraction. RSA improves by 4.2% (up to 5.5% for 512-bit RSA). The GMPbench results expose a structural property of multi-precision libraries: speeding up the primitives is not an isolated issue. A focused improvement at the bottom of the stack propagates broadly across the suite. OpenSSL. Figure 5 shows DoTSSL’s throughput improvements on OpenSSL’s speed benchmark for RSA, RSA KEM, FFDH, and DSA. DoTSSL consistently improves throughput across all tested key sizes and operations. For RSA (1024– 7680-bit keys), improvements average 3.2% across sign, verify, encrypt, and decrypt. FFDH shows the most consistent gains, averaging 4.4% across group sizes (up to 5.9% for 4096-bit groups). DSA improves by up to 5.2% (verify, 2048-bit). DoT’s Contribution. To understand the contribution of DoT’s addition, subtraction, and multiplication routines to the observed application-level speedups, we performed an analysis of the cycles spent (%) in the GMPbench and OpenSSL speed benchmark workloads for DoT’s replaced routines:

Conclusion

Large-number addition and subtraction remain unaccelerated by SIMD in every major multi-precision library, despite being the foundation of cryptography and high-precision computing. We present DigitsOnTurbo (DoT): a 4-phase SIMD add/sub algorithm that breaks the sequential carry chain by deferring rare cascades to a slow path, and a “vertical and crosswise”-based multiplication routine that exposes all 𝑚 2 partial products as independent IFMA work. Both are implemented in portable C with AVX-512 intrinsics, no hand-tuned assembly. On Emerald Rapids (𝑤=8), DoT add/sub reach 1.85×/1.84× over scalar add-with-carry and DoT multiplication reaches 2.3× over the AVX-512IFMA 10

Leveraging SIMD for Accelerating Large-number Arithmetic <1.5% +2.1%

128 512 8K 128K 2048K 128x128 512x512 Multiply 8Kx8K 128Kx128K 2048Kx2048K 15Kx10K 20Kx10K 30Kx10K 16384Kx512 16384Kx256K

<1.5% <1.5% <1.5% +2.6%

128 512

+10.3% +17.5%

GCD +28.8%

8K 128K

<1.5% +35.6%

<1.5% <1.5%

128

+48.7% +23.1% +18.6% +16.0%

GCDext

512

-4.1%

GCDext

+3.1%

GCD

+3.1% +8.4%

Divide

+6.1%

128K

+4.0%

+4.2%

RSA

+13.1%

1024K

+2.8% +9.9%

+13.3%

+15.3%

Multiply

+13.7%

1024K

+22.9% <1.5% <1.5% <1.5% <1.5%

8K÷32 8K÷64 8K÷128 8K÷4K Divide128K÷64K 8192K÷4096K 8K÷8064 16384K÷256K

RSA +17.7%

Pi 100K

+12.8%

−5

(a) Mul/Div Improvement (%)

+7.3%

total

+7.8%

+19.3%

1000K

+18.6%

+8.7%

app base

+8.0%

10K

+36.4%

1K 2K

<1.5%

−5

+5.5% +1.7% +5.5%

512

−5

(b) GCD/RSA/Pi Improvement (%)

Figure 4. DoTMP’s score (throughput) improvement over GMP in GMPbench. Sign/s Verify/s

Improvement (%)

Encrypt/Encaps Decrypt/Decaps

Keygen (op/s) 6

0 1024

2048

3072

4096

7680

0 1024

Key Size (bits) (a) RSA

2048

3072

4096

7680

Sign/s Verify/s

0 2048

Key Size (bits) (b) RSA KEM

3072

4096

6144

8192

1024

Group Size (bits) (c) FFDH

2048

Key Size (bits) (d) DSA

Figure 5. DoTSSL throughput improvement (%) over OpenSSL for RSA (sign/verify/encrypt/decrypt), RSA KEM (encaps/decaps), FFDH (keygen), and DSA (sign/verify) across standard key and group sizes. 512

GCD

512×512

1024

128K

2048

RSA 3072

4096

8K×8K 15K×10K

GCDext

20K×10K

7680

128K 1M

Multiply30K×10K 128K 128K×128K

512

RSA

2M×2M

DSA

1024 2048

16M×512

dot_mul_4x4 dot_add_words dot_sub_words

16M×256K 128K÷64K

10K

2048

Divide 8M÷4M

Pi 100K

FFDH 3072

16M÷256K

4096

1M 0

Cycles Spent (%) in DoT Routines

(a) GMPbench (Mul, Div)

(b) GMPbench (GCD, RSA, Pi)

Figure 6. Cycle spent (%) by DoT’s dot_add_words, dot_sub_words, and dot_mul_4x4 routines in GMPbench and OpenSSL speed workloads, measured via perf. We omitted a handful of cases in the GMPbench (e.g., lower sized mul, div and gcd) since they spend zero cycles in DoT routines. baseline for 256-bit operands. Integrated into GMP (DoTMP) and OpenSSL (DoTSSL), these gains propagate end-to-end: GMPbench’s overall score improves by 7.8% (𝜋 up to +19.3%,

multiply up to +48.7%), and OpenSSL throughput for RSA, FFDH, and DSA improves by up to 5.9%. The broader takeaway is that SIMD adoption in large-number arithmetic has 11

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

been blocked by algorithmic structure, not by hardware. Restructuring around independent and data-parallel operations is what unlocks the gains offered by SIMD.

[21] Laurent-Stéphane Didier, Nadia El Mrabet, Léa Glandus, and JeanMarc Robert. 2024. Truncated multiplication and batch software SIMD AVX512 implementation for faster Montgomery multiplications and modular exponentiation. IACR Communications in Cryptology 1, 3 (2024). doi:10.62056/a3txl86bm [22] Whitfield Diffie and Martin E. Hellman. 2022. New Directions in Cryptography (1 ed.). Association for Computing Machinery, New York, NY, USA, 365–390. https://doi.org/10.1145/3549993.3550007 [23] Mozilla JS Docs. 2025. BigInt - JavaScript | MDN — developer.mozilla.org. https://developer.mozilla.org/en-US/docs/Web/ JavaScript/Reference/Global_Objects/BigInt. [Accessed 12-03-2025]. [24] Takuya Edamatsu and Daisuke Takahashi. 2018. Acceleration of Large Integer Multiplication with Intel AVX-512 Instructions. In 2018 IEEE 20th International Conference on High Performance Computing and Communications; IEEE 16th International Conference on Smart City; IEEE 4th International Conference on Data Science and Systems (HPCC/SmartCity/DSS). 211–218. doi:10.1109/HPCC/SmartCity/DSS. 2018.00059 [25] Takuya Edamatsu and Daisuke Takahashi. 2019. Accelerating Large Integer Multiplication Using Intel AVX-512IFMA. In Algorithms and Architectures for Parallel Processing: 19th International Conference, ICA3PP 2019, Melbourne, VIC, Australia, December 9–11, 2019, Proceedings, Part I (Melbourne, VIC, Australia). Springer-Verlag, Berlin, Heidelberg, 60–74. doi:10.1007/978-3-030-38991-8_5 [26] Takuya Edamatsu and Daisuke Takahashi. 2023. Efficient Large Integer Multiplication with Arm SVE Instructions. In Proceedings of the International Conference on High Performance Computing in Asia-Pacific Region (Singapore, Singapore) (HPCAsia ’23). Association for Computing Machinery, New York, NY, USA, 9–17. doi:10.1145/3578178.3578193 [27] Andres Erbsen, Jade Philipoom, Jason Gross, Robert Sloan, and Adam Chlipala. 2020. Simple High-Level Code For Cryptographic Arithmetic: With Proofs, Without Compromises. SIGOPS Oper. Syst. Rev. 54, 1 (Aug. 2020), 23–30. doi:10.1145/3421473.3421477 [28] FLINT Development Team. 2025. FLINT: Fast Library for Number Theory — flintlib.org. https://flintlib.org/. [Accessed 05-05-2025]. [29] M.J. Flynn. 1966. Very high-speed computing systems. Proc. IEEE 54, 12 (1966), 1901–1909. doi:10.1109/PROC.1966.5273 [30] Agner Fog. 2025. 4. Instruction tables Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD, and VIA CPUs. https://www.agner.org/optimize/instruction_tables.pdf. [Accessed 14-09-2025]. [31] Gerhard Frey. 2010. The arithmetic behind cryptography. Notices of the AMS 57, 3 (2010), 366–374. [32] GCC 2025. GCC, the GNU Compiler Collection - GNU Project — gcc.gnu.org. https://gcc.gnu.org/. [Accessed 24-03-2026]. [33] GMPbench. 2025. GMPbench results — gmplib.org. https://gmplib. org/gmpbench. [Accessed 21-03-2025]. [34] GNU Project. 1991. The GNU MP Bignum Library — gmplib.org. https://gmplib.org/. [Accessed 03-03-2025]. [35] Shay Gueron and Vlad Krasnov. 2012. Software Implementation of Modular Exponentiation, Using Advanced Vector Instructions Architectures. In Arithmetic of Finite Fields, Ferruh Özbudak and Francisco Rodríguez-Henríquez (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 119–135. [36] Shay Gueron and Vlad Krasnov. 2015. Fast prime field elliptic-curve cryptography with 256-bit primes. Journal of Cryptographic Engineering 5, 2 (2015), 141–151. doi:10.1007/s13389-014-0090-x [37] Shay Gueron and Vlad Krasnov. 2016. Accelerating Big Integer Arithmetic Using Intel IFMA Extensions. In 2016 IEEE 23nd Symposium on Computer Arithmetic (ARITH). 32–38. doi:10.1109/ARITH.2016.22 [38] Martin E. Hellman. 1979. The Mathematics of Public-Key Cryptography. Scientific American 241, 2 (1979), 146–157. http://www.jstor.org/ stable/24965269

References [1] 2017. Intel Advanced Vector Extensions 512 (Intel AVX-512) Overview — intel.com. https://www.intel.com/content/www/us/en/architectureand-technology/avx-512-overview.html. [Accessed 16-09-2025]. [2] 2026. Simd Library — ermig1979.github.io. https://ermig1979.github. io/{S}imd/. [Accessed 03-04-2026]. [3] Advanced Micro Devices, Inc. 2025. Leadership HPC Performance with 5th Generation AMD EPYC Processors. https://www.amd.com/en/blogs/2025/leadership-hpc-performancewith-5th-generation-amd.html. [4] Arm ADC 2022. Documentation; Arm Developer — developer.arm.com. https://developer.arm.com/documentation/ddi0602/ 2022-06/Base-Instructions/ADC--Add-with-Carry-. [Accessed 18-092025]. [5] Arm PL 2025. Arm Performance Libraries — developer.arm.com. https://developer.arm.com/{T}ools%20and%20{S}oftware/{A}rm% 20{P}erformance%20{L}ibraries. [Accessed 25-03-2026]. [6] Arm SVE2 2022. Documentation; Arm Developer — developer.arm.com. https://developer.arm.com/documentation/102340/ latest/SVE2-architecture-fundamentals. [Accessed 19-09-2025]. [7] D.H. Bailey. 2005. High-precision floating-point arithmetic in scientific computation. Computing in Science & Engineering 7, 3 (2005), 54–61. doi:10.1109/MCSE.2005.52 [8] D.H. Bailey, R. Barrio, and J.M. Borwein. 2012. High-precision computation: Mathematical physics and dynamics. Appl. Math. Comput. 218, 20 (2012), 10106–10121. doi:10.1016/j.amc.2012.03.087 [9] David H. Bailey and Jonathan M. Borwein. 2015. High-Precision Arithmetic in Mathematical Physics. Mathematics 3, 2 (2015), 337–367. doi:10.3390/math3020337 [10] Elaine Barker. 2020. Recommendation for Key Management: Part 1 – General. https://doi.org/10.6028/NIST.SP.800-57pt1r5. [Accessed 13-03-2025]. [11] O. J. Bedrij. 1962. Carry-Select Adder. IRE Transactions on Electronic Computers EC-11, 3 (1962), 340–346. doi:10.1109/IRETELC.1962. 5407919 [12] Clifton Haider Benjamin Buhrow, Barry Gilbert. 2021. Parallel modular multiplication using 512-bit advanced vector instructions - Journal of Cryptographic Engineering — link.springer.com. https://link. springer.com/article/10.1007/s13389-021-00256-9. doi:10.1007/s13389021-00256-9 [Accessed 08-09-2025]. [13] Andrew D Booth. 1951. A signed binary multiplication technique. The Quarterly Journal of Mechanics and Applied Mathematics 4, 2 (1951), 236–240. [14] Brent and Kung. 1982. A regular layout for parallel adders. IEEE transactions on Computers 100, 3 (1982), 260–264. [15] Richard Brent and Paul Zimmermann. 2010. Modern Computer Arithmetic. Cambridge University Press, USA. [16] Lin Chao. 1999. Intel Technology Journal Q2. https://www.intel.com/ content/dam/www/public/us/en/documents/research/1999-vol03iss-2-intel-technology-journal.pdf. [Accessed 16-03-2025]. [17] Neil Coffey. 2025. RSA key lengths — javamex.com. https://www. javamex.com/tutorials/cryptography/rsa_key_length.shtml. [Accessed 12-03-2025]. [18] P. G. Comba. 1990. Exponentiation cryptosystems on the IBM PC. IBM Systems Journal 29, 4 (1990), 526–538. doi:10.1147/sj.294.0526 [19] Intel Cooperation. 2000. Using Streaming SIMD Extensions (SSE2) to Perform Big Multiplications. Technical Report. Technical Report. [20] Luigi Dadda. 1965. Some schemes for parallel multipliers. Alta frequenza 34 (1965), 349–356. 12

Leveraging SIMD for Accelerating Large-number Arithmetic [39] John L. Hennessy and David A. Patterson. 2012. Computer Architecture: A Quantitative Approach (5th ed.). Morgan Kaufmann / Elsevier. [40] Mike Housch. 2025. The Current Encryption Landscape: The Need For 3072-Bit Keys — forbes.com. https://www.forbes.com/councils/ forbestechcouncil/2024/02/23/the-current-encryption-landscapethe-need-for-3072-bit-keys/. [Accessed 12-03-2025]. [41] The MathWorks Inc. 2022. Symbolic Math Toolbox. https://in. mathworks.com/products/symbolic.html [42] Intel. 2025. Intel® Advanced Vector Extensions 10.1 (Intel® AVX10.1) Architecture Specification — intel.com. https://www.intel.com/ content/www/us/en/content-details/848455/intel-advanced-vectorextensions-10-1-intel-avx10-1-architecture-specification.html. [Accessed 30-04-2025]. [43] Intel AVX2 2021. Intel; Advanced Vector Extensions 2 (Intel AVX-2) - 009 - ID:655258; Processors — edc.intel.com. https://edc.intel.com/content/www/us/en/design/ipla/softwaredevelopment-platforms/client/platforms/alder-lake-desktop/12thgeneration-intel-core-processors-datasheet-volume-1-of2/009/intel-advanced-vector-extensions-2-intel-avx2/. [Accessed 16-09-2025]. [44] Intel Corporation 2024. Intel® 64 and IA-32 Architectures Optimization Reference Manual. Intel Corporation. Volume 1, Document 248966050, April 2024. See Chapter 18 (Software Optimization for Intel AVX512 Instructions) for general pipeline, dependency, and accumulator guidance on fused-multiply-accumulate style operations; Chapter 21.4 (or Chapter 19.4 in some printings) gives the functional description of AVX512_IFMA (VPMADD52) instructions and their intended use in big-number multiplication.. [45] Intel MKL 2025. Accelerate Fast Math with Intel® oneAPI Math Kernel Library — intel.com. https://www.intel.com/content/www/us/ en/developer/tools/oneapi/onemkl.html. [Accessed 25-03-2026]. [46] Intel SDM 2025. Manuals for Intel® 64 and IA-32 Architectures — intel.com. https://www.intel.com/content/www/us/en/developer/ articles/technical/intel-sdm.html. [Accessed 18-09-2025]. [47] IntelIntrins. 2024. Intel® Intrinsics Guide — intel.com. https://www. intel.com/content/www/us/en/docs/intrinsics-guide/index.html. [Accessed 05-03-2025]. [48] Fredrik Johansson. 2025. mpmath - Python library for arbitraryprecision floating-point arithmetic — mpmath.org. https://mpmath. org/. [Accessed 12-03-2025]. [49] Don Johnson, Alfred Menezes, and Scott Vanstone. 2001. The Elliptic Curve Digital Signature Algorithm (ECDSA). Int. J. Inf. Secur. 1, 1 (Aug. 2001), 36–63. doi:10.1007/s102070100002 [50] Anatolii Karatsuba. 1963. Multiplication of multidigit numbers on automata. In Soviet physics doklady, Vol. 7. 595–596. [51] Anastasis Keliris and Michail Maniatakos. 2014. Investigating large integer arithmetic on Intel Xeon Phi SIMD extensions. In 2014 9th IEEE International Conference on Design & Technology of Integrated Systems in Nanoscale Era (DTIS). 1–6. doi:10.1109/DTIS.2014.6850661 [52] Donald E Knuth. 1997. The Art of Computer Programming, Volume 2: Seminumerical Algorithms (third ed.). Addison-Wesley Professional, Boston. [53] Peter M. Kogge and Harold S. Stone. 1973. A Parallel Algorithm for the Efficient Solution of a General Class of Recurrence Equations. IEEE Trans. Comput. 22, 8 (Aug. 1973), 786–793. doi:10.1109/TC.1973. 5009159 [54] Feng Liu, Qingping Tan, and Gang Chen. 2010. Formal proof of prefix adders. Mathematical and Computer Modelling 52, 1 (2010), 191–199. doi:10.1016/j.mcm.2010.02.008 [55] LLVM Overflow 2025. LLVM Language Reference Manual; LLVM 22.0.0git documentation — llvm.org. https://llvm.org/docs/LangRef. html. [Accessed 18-09-2025]. [56] O. L. Macsorley. 1961. High-Speed Arithmetic in Binary Computers. Proceedings of the IRE 49, 1 (1961), 67–91. doi:10.1109/JRPROC.1961.

287779 Vedic Mathematics. [57] Bharati Krsna Tirthji Maharaj. 1992. https://archive.org/details/vedic-mathematics-bharati-krishnatirth-ji-maharaj/page/n7/mode/2up. [Accessed 05-03-2025]. [58] Linux man pages. 2024. perf_event_open(2) - Linux manual page — man7.org. https://www.man7.org/linux/man-pages/man2/perf_ event_open.2.html. [Accessed 20-03-2025]. [59] Makoto Matsumoto and Takuji Nishimura. 1998. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Trans. Model. Comput. Simul. 8, 1 (Jan. 1998), 3–30. doi:10.1145/272991.272995 [60] Maxima. 2025. Maxima – GPL CAS based on DOE-MACSYMA — maxima.sourceforge.io. https://maxima.sourceforge.io/. [Accessed 12-03-2025]. [61] Victor S. Miller. 1986. Use of Elliptic Curves in Cryptography. In Advances in Cryptology — CRYPTO ’85 Proceedings, Hugh C. Williams (Ed.). Springer Berlin Heidelberg, Berlin, Heidelberg, 417–426. doi:10. 1007/3-540-39799-X_31 [62] Mala Saraswathy Nataraj and Michael O. J. Thomas. 2006. Expansion of binomials and factorisation of quadratic expressions: Exploring a Vedic method. Australian Senior Mathematics Journal 20, 2 (2006), 8–17. [63] Linux on IBM Systems. 2025. Common Cryptographic Architecture (CCA): ECC key token — ibm.com. https://www.ibm.com/docs/en/ linux-on-systems?topic=formats-ecc-key-token. [Accessed 13-032025]. [64] OpenBLAS 2025. OpenBLAS : An optimized BLAS library — openmathlib.org. http://www.openmathlib.org/{O}pen{B}{L}{A}{S}. [Accessed 25-03-2026]. [65] OpenSSL rsaz 2025. Openssl RSAZ. https://github.com/openssl/ openssl/blob/master/crypto/bn/rsaz_exp_x2.c. [Accessed 06-09-2025]. [66] OpenSSL Software Foundation. 2025. OpenSSL — openssl.org. https: //www.openssl.org/. [Accessed 05-05-2025]. [67] Gabriele Paoloni. 2010. How to Benchmark Code Execution Times on Intel IA-32 and IA-64 Instruction Set Architectures. Intel White Paper. [Accessed 21-03-2025]. [68] GNU Project. 2025. The GNU MPFR Library — mpfr.org. https://www. mpfr.org/. [Accessed 12-03-2025]. [69] Pengchang Ren, Reiji Suda, and Vorapong Suppakitpaisarn. 2023. Efficient Additions and Montgomery Reductions of Large Integers for SIMD. In 2023 IEEE 30th Symposium on Computer Arithmetic (ARITH). 48–59. doi:10.1109/ARITH58626.2023.00034 [70] R. L. Rivest, A. Shamir, and L. Adleman. 1978. A method for obtaining digital signatures and public-key cryptosystems. Commun. ACM 21, 2 (Feb. 1978), 120–126. doi:10.1145/359340.359342 [71] SageMath. 2025. SageMath Mathematical Software System - Sage — sagemath.org. https://www.sagemath.org/. [Accessed 12-03-2025]. [72] Arnold Schönhage and Volker Strassen. 1971. Fast multiplication of large numbers. Computing 7 (1971), 281–292. [73] GNU MP SIMD. 2025. Assembly SIMD Instructions (GNU MP 6.3.0) — gmplib.org. https://gmplib.org/manual/Assembly-SIMD-Instructions. [Accessed 12-03-2025]. [74] J. Sklansky. 1960. Conditional-Sum Addition Logic. IRE Transactions on Electronic Computers EC-9, 2 (1960), 226–231. doi:10.1109/TEC.1960. 5219822 [75] SSL Support Team. 2025. New Minimum RSA Key Size for Code Signing Certificates - SSL.com — ssl.com. https://www.ssl.com/blogs/newminimum-rsa-key-size-for-code-signing-certificates/. [Accessed 13-03-2025]. [76] Mikko Tommila. 2025. Apfloat - Arbitrary precision library for Java and C++, applets and calculator. http://www.apfloat.org/. [Accessed 12-03-2025]. [77] Andrei L Toom. 1963. The complexity of a scheme of functional elements realizing the multiplication of integers, published in Soviet 13

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel Math (translations of Dokl. Adad. Nauk. SSSR), 4. [78] Daniel Towner. 2022. Intel Advanced Vector Extensions 512 (Intel AVX-512) - Permuting Data Within and Between AVX Registers. https://builders.intel.com/docs/networkbuilders/intel-avx-512permuting-data-within-and-between-avx-registers-technologyguide-1668169807.pdf. [Accessed 16-03-2025]. [79] Christopher S Wallace. 2006. A suggestion for a fast multiplier. IEEE Transactions on electronic Computers 1 (2006), 14–17. [80] Lynn West. 2011. An Introduction to Various Multiplication Strategies. https://www.educator.com/classroom/users/h/highgater/961_ Many_Ways_to_Multiply.pdf. [81] Wolfram. 2025. Wolfram Mathematica: Modern Technical Computing — wolfram.com. https://www.wolfram.com/mathematica/. [Accessed 12-03-2025]. [82] Alexander Yee. 2019. Integer Addition and Carryout — numberworld.org. http://numberworld.org/y-cruncher/internals/addition. html. [Accessed 03-03-2025]. [83] Alexander J. Yee. 2025. y-cruncher - A Multi-Threaded Pi Program — numberworld.org. http://numberworld.org/y-cruncher/. [Accessed 06-05-2025].

Leveraging SIMD for Accelerating Large-number Arithmetic

Hardware and Library Context for Limb-level Arithmetic

Algorithm 4: Karatsuba Multiplication (or, ToomCook 2-way) Data: Integers 𝐴, 𝐵 with 𝑚 limbs Result: Product 𝑃 = 𝐴 · 𝐵 1 Function: MPN-KARATSUBA(𝐴, 𝐵, 𝑚) 2 begin 3 if 𝑚 ≤ 𝜃 then 4 return MPN-MUL-BASECASE(𝐴, 𝐵, 𝑚); 5 end 6 𝑘 ← 𝑚/2, 𝑏 ← 2𝑘 ·bits_per_limb ; 7 𝐴𝐻 , 𝐴𝐿 ← split 𝐴 at 𝑘; 𝐵𝐻 , 𝐵𝐿 ← split 𝐵 at 𝑘; 8 𝑃 1 ← MPN-KARATSUBA(𝐴𝐻 , 𝐵𝐻 , 𝑘 ); 9 𝑃 0 ← MPN-KARATSUBA(𝐴𝐿 , 𝐵𝐿 , 𝑘 ); 10 𝑠𝑥 , 𝑠 𝑦 ← sign of(𝐴𝐻 − 𝐴𝐿 ), sign of(𝐵𝐻 − 𝐵𝐿 ); 11 𝑃 Δ ← MPN-KARATSUBA( |𝐴𝐻 − 𝐴𝐿 |, |𝐵𝐻 − 𝐵𝐿 |, 𝑘 ); 12 𝑃 2 ← (𝑠𝑥 · 𝑠 𝑦 )𝑃 Δ ; 13 return (𝑏 2 + 𝑏 )𝑃1 − 𝑏𝑃2 + (𝑏 + 1)𝑃 0 ; 14 end

This appendix expands on the scalar hardware support and library-level context that Section 2.2 only summarizes. ALU width and word-local optimizations. On most CPUs the arithmetic logic unit (ALU) handles fixed-size operands of 32 or 64 bits. Within a single word, hardware designs such as carry-lookahead [56] and carry-select logic [11], as well as parallel-prefix trees (Kogge-Stone [53], Sklansky [74], Brent-Kung [14]), reduce the bit-level carry delay from ripple-carry 𝑂 (𝑛) toward prefix-style 𝑂 (log 𝑛) [39]. Modern integer pipelines hide much of this work in practice. The catch is that all of these optimizations apply only within a single machine word, not across the multiple words that make up a large-number operand. Limb-by-limb addition with hardware carry chains. For large-numbers, the operation has to be decomposed into a sequence of word-sized limb additions [15]. For example, adding two 512-bit numbers on a 64-bit architecture takes eight 64-bit limb additions. Each one is fast on its own, but a sequential dependency creeps in: the carry-out of each addition has to flow into the next as carry-in, which brings the linear 𝑂 (𝑚) delay back at the limb level. The textbook software approach computes 𝑆𝑖 = 𝐴𝑖 +𝐵𝑖 +𝐶𝑖 −1 [15, 52]; Algorithm 3 (MPN-ADD-M) shows GMP’s realization of this idea. To avoid manual carry detection and propagation, low-level assembly uses the add-with-carry chaining technique [4, 46]: the least significant limb addition starts with the ADD instruction (which sets the hardware carry flag), and subsequent limbs use ADC instructions that add both the limb operands and the incoming carry flag. Subtraction uses SBB instructions in the same way. This pattern is widely deployed in hand-optimized GMP and OpenSSL assembly, because mainstream compilers still do not reliably synthesize equivalent carry chains from straightforward high-level C loops [32, 55].

Hardware multipliers and base-case multiplication routines. Arbitrary-precision multiplication builds on the same fixed-size hardware multipliers. A modern 64×64-bit unit uses partial-product generation, booth’s encoding [13], compressor-tree reduction (typically Wallace or Dadda trees [20, 79]), and a final carry-propagate addition, all within a few cycles of latency [30, 39]. Many implementations also use Booth-style recoding to cut down on partial products before reduction [39]. These datapath optimizations help single-limb multiply throughput but do nothing about the cross-limb accumulation dependencies that arise in arbitraryprecision multiplication. For small sizes (𝑚 = 2, . . . , 8), libraries primarily use schoolbook and Comba [15, 18, 52]: schoolbook accumulates row-wise, while Comba accumulates column-wise to shorten carry chains. Assembly implementations commonly use mulx, which computes the full 128-bit product of two 64-bit limbs and returns both halves in a single instruction.

Algorithm 3: Sequential Addition (GMP-style)

Recursive multiplication thresholds. For larger operands, most libraries fall back to Karatsuba’s [15, 50] divide-andconquer approach (Algorithm 4), which uses three half-size multiplications instead of four and brings the complexity down to 𝑂 (𝑚 log2 3 ) ≈ 𝑂 (𝑚 1.585 ). GMP further switches to Toom-Cook [77] (a generalization of Karatsuba) and to FFTbased multiplication [72] for very large operands, while OpenSSL typically only goes as far as Karatsuba beyond the base case. The exact switch points depend on the implementation, the operand size, and the target architecture. When the recursion reaches its base case (Line 3), the library invokes schoolbook or Comba. Crucially, because Karatsuba and Toom-Cook trade multiplications for additions and subtractions at every recursive node, the aggregate add/sub work grows quickly with recursion depth. So a faster base-case

Data: Integers 𝐴, 𝐵 each with 𝑚 limbs Result: Sum 𝑅 and final carry 𝑐𝑜𝑢𝑡 1 // Sequential limb addition with carry propagation 2 Function: MPN-ADD-M(𝑅, 𝐴, 𝐵, 𝑚) 3 begin 4 𝑐𝑖𝑛 ← 0; 5 for 𝑖 ← 0 to 𝑚 − 1 step 1 do 6 𝑅 [𝑖 ] ← 𝐴[𝑖 ] + 𝑐𝑖𝑛 ; 7 𝑐𝑜𝑢𝑡 ← (𝑅 [𝑖 ] < 𝑐𝑖𝑛 ) ? 1 : 0; 8 𝑅 [𝑖 ] ← 𝑅 [𝑖 ] + 𝐵 [𝑖 ]; 9 𝑐𝑜𝑢𝑡 ← 𝑐𝑜𝑢𝑡 + ( (𝑅 [𝑖 ] < 𝐵 [𝑖 ] ) ? 1 : 0); 10 end 11 return 𝑐𝑜𝑢𝑡 ; 12 end

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

multiply and a faster add/sub both feed into every higherlevel operation built on top: division, modular exponentiation, and computations like 𝜋.

To show that if 𝑠 < 𝑎 and 𝑠 < 𝑏, then a carry is generated: Assume 𝑠 < 𝑎 and no carry is generated. As no carry is generated, 0 ≤ 𝑠 = 𝑎 + 𝑏 < 2𝑛 . Since 𝑏 ≥ 0, 𝑠 ≥ 𝑎, which is a contradiction. Similarly, if 𝑠 < 𝑏, a carry has to be generated. □

Proof of Correctness and Carry Probability Analysis

Theorem B.2 (Correctness of DoT-Addition). Suppose two large integers 𝐴 and 𝐵,

This appendix formalizes the correctness of DoT’s addition and multiplication. Also, it analyzes the probability of carry generation in large-number addition, which motivates DoT’s design choice to isolate carry handling into a separate phase that only triggers on rare carry cascades.

𝐴=

𝑛−1 ∑︁

𝐵=

𝑛−1 ∑︁

𝐵 𝑗 𝑋 𝑗 = 𝐵𝑛−1𝑋 𝑛−1 + · · · + 𝐵 1𝑋 1 + 𝐵 0𝑋 0

𝑗=0

where 𝑋 = 2𝑘 (where 𝑘 is the bit size of a limb), the number of limbs for the two integers is 𝑛, and 0 ≤ 𝐴𝑖 , 𝐵 𝑗 < 𝑋 . Then, if DoT-Add-Words(𝐴, 𝐵, 𝑛) = (𝑆, 𝑐 out ), then

𝐴𝑖 𝑋 𝑖 = 𝐴𝑛−1𝑋 𝑛−1 + · · · + 𝐴1𝑋 1 + 𝐴0

𝑆=

𝑖=0

𝐵=

𝐴𝑖 𝑋 𝑖 = 𝐴𝑛−1𝑋 𝑛−1 + · · · + 𝐴1𝑋 1 + 𝐴0𝑋 0

𝑖=0

Proof of correctness of DoT’s addition algorithm. Let two large integers 𝐴 and 𝐵 be represented in radix 𝑋 = 2𝑘 (where 𝑘 is the bit size of a limb) and suppose that the number of limbs for the two integers is 𝑛: 𝐴=

𝑛−1 ∑︁

(𝑟𝑖 𝑋 𝑖 )

𝑖=0

𝐵 𝑗 𝑋 𝑗 = 𝐵𝑛−1𝑋 𝑛−1 + · · · + 𝐵 1𝑋 1 + 𝐵 0

such that 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋, and 𝑐 out = 𝑐𝑛

𝑗=0

where 0 ≤ 𝐴𝑖 , 𝐵 𝑗 < 𝑋 . The sum of these two numbers using schoolbook addition is given by: 𝑆=

𝑛−1 ∑︁

Proof. The proof proceeds by induction on the DoT-AddWords algorithm. Base case (for n = 1 limb): For the base case, this reduces to the addition of two 64-bit integers. In Phase 1, we compute 𝑟 0 = 𝐴0 + 𝐵 0 . In Phase 2, if 𝑟 0 < 𝐴0 , then 𝑐 = 1 else 𝑐 = 0. As 𝑤 = 1, 𝑐𝑜𝑢𝑡 = 𝑐 ≫ 0, i.e., 𝑐𝑜𝑢𝑡 = 𝑐. In Phase 3, 𝑟 0 = 𝐴0 + 𝐵 0 +𝑐 0 . As 𝑐 0 = 0, 𝑟 0 = 𝐴0 + 𝐵 0 . Because, 𝑟 0 is unchanged, Phase 4 is not triggered. Thus, the final value is 𝑆 = 𝐴0 + 𝐵 0 and 𝑐𝑜𝑢𝑡 = 𝑐. From Lemma B.1, we know that if carry is generated 𝑆 < 𝐴0 and 𝑆 < 𝐵 0 . Inductive case (w.l.o.g, we show for 0 < 𝑛 < 8). IH: for 𝑛 limbs, 𝑛−1 ∑︁ 𝑆′ = (𝑟𝑖 𝑋 𝑖 )

(𝑠𝑖 𝑋 𝑖 ) + 𝑐𝑛 𝑋 𝑛

𝑖=0

such that 𝑠𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋ To show that Algorithm 1 is correct, we need to show that the output of the algorithm is same as the sum computed above or 𝑛−1 ∑︁ 𝑆= (𝑠𝑖 𝑋 𝑖 ) 𝑖=0

and 𝑐𝑜𝑢𝑡 = 𝑐𝑛 . Since the maximum value of a limb (𝐴𝑖 , 𝐵𝑖 ) is 𝑋 − 1, the maximum sum of two limbs and a carry is: (𝑋 − 1) + (𝑋 − 1) + 1 = 2𝑋 − 1. This means that 𝑐𝑖+1 can never exceed ⌊(2𝑋 − 1)/𝑋 ⌋ = 1. Also, when 𝑠𝑖 < 𝐴𝑖 or 𝑠𝑖 < 𝐵𝑖 , it means that a carry is generated because the sum of two numbers is always greater than the individual numbers.

𝑖=0

such that 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋, and 𝑐 out = 𝑐𝑛 . To show: for (𝑛 + 1)𝑡ℎ limb, 𝑛 ∑︁ 𝑆= (𝑟𝑖 𝑋 𝑖 )

Lemma B.1 (Carry generation exceeds limb size). The sum of two 𝑛-bit integers 𝑎 and 𝑏 generates a carry, if and only if the least significant 𝑛 bits of the sum are less than both 𝑎 and 𝑏. In other words, if 𝑎 + 𝑏 ≥ 𝑎, then 𝑎 + 𝑏 ≥ 𝑏 and no carry-bit is generated.

𝑖=0

such that 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋, and 𝑐 out = 𝑐𝑛+1 .

Proof. Let 𝑎 and 𝑏 be 𝑛-bit integers such that 0 ≤ 𝑎, 𝑏 < 2𝑛 , and 𝑠 = 𝑎 + 𝑏. A carry occurs if 𝑠 ≥ 2𝑛 . Since 𝑎, 𝑏 < 2𝑛 , the maximum value of 𝑠 is 2𝑛 − 1 + 2𝑛 − 1 = 2𝑛+1 − 2, which is strictly less than 2𝑛+1 . As the sum is 2𝑛 ≤ 𝑠 < 2𝑛+1 , the last 𝑛-bits of the sum are given by 𝑠 mod 2𝑛 , which is 𝑠 − 2𝑛 . Since 𝑎 and 𝑏 are less than 2𝑛 , the sum is less than both 𝑎 and 𝑏.

𝑆 = (𝑟𝑛 𝑋 𝑛 ) +

𝑛−1 ∑︁

(𝑟𝑖 𝑋 𝑖 )

𝑖=0

𝑆 = (𝑟𝑛 𝑋 𝑛 ) + 𝑆 ′ such that 𝑐𝑛 = ⌊(𝐴𝑛−1 + 𝐵𝑛−1 + 𝑐𝑛−1 )/𝑋 ⌋ From Phase 1, we have ∀ (0 ≤ 𝑖 ≤ 𝑛). 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 ) 16

Leveraging SIMD for Accelerating Large-number Arithmetic

In Phase 2, if 𝑟𝑛 < 𝐴𝑛 , then 𝑐𝑛 = 1 else 𝑐𝑛 = 0. As 𝑤 = 𝑛, the MSBs. 𝑐𝑜𝑢𝑡 = 𝑐 ≫ 𝑛, i.e., 𝑐𝑜𝑢𝑡 = 𝑐𝑛 . From IH, we have ∀ (0 ≤ 𝑛−1 𝑛−1 ∑︁ ∑︁ 𝑖 < 𝑛). 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = 𝐴×𝐵 = 𝐴𝑖 𝑋 𝑖 × 𝐵𝑗𝑋 𝑗 ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋. Thus, the new 𝑐𝑛 = ⌊(𝐴𝑛−1 + 𝐵𝑛−1 + 𝑐𝑛−1 )/𝑋 ⌋, 𝑖=0 𝑗=0 which is 1 if 𝑟𝑛−1 < 𝐴𝑛−1 else 0. =(𝐴𝑛−1𝑋 𝑛−1 + · · · + 𝐴1𝑋 1 + 𝐴0𝑋 0 )× In Phase 3, we have 𝑟𝑛 = 𝐴𝑛 + 𝐵𝑛 + 𝑐𝑛 . If this is more than 𝑋 , 𝐵𝑛−1𝑋 𝑛−1 + · · · + 𝐵 1𝑋 1 + 𝐵 0𝑋 0 we get 𝑟𝑛 = (𝐴𝑛 + 𝐵𝑛 +𝑐𝑛 ) (mod 𝑋 ), which means a carry is further generated, i.e., (𝐴𝑛 + 𝐵𝑛 + 𝑐𝑛 ) (mod 𝑥) < (𝐴𝑛 + 𝐵𝑛 ) =(𝐴𝑛−1𝑋 𝑛−1 .𝐵𝑛−1𝑋 𝑛−1 + 𝐴𝑛−1𝑋 𝑛−1 .𝐵𝑛−2𝑋 𝑛−2 + ...+ or ⌊(𝐴𝑛 + 𝐵𝑛 + 𝑐𝑛 )/𝑋 ⌋ = 1 (from Lemma B.1). The Phase 4 𝐴0𝑋 0 .𝐵 0𝑋 0 triggers when there is some carry generated in one of the 𝑛−1 ∑︁ 𝑛−1 ∑︁ limbs. The correctness of Phase 4 and the final carry genera= 𝐴𝑖 .𝐵 𝑗 𝑋 𝑖+𝑗 tion follows from the correctness of Kogge-Stone Adder [54]. 𝑖=0 𝑗=0 Thus, we have for 𝑛 + 1 limbs, (1) 𝑆 = 𝑟𝑛 𝑋 𝑛 +

𝑛−1 ∑︁

DoT-Mul-Words(𝐴, 𝐵, 𝑛, 𝑘)

(𝑟𝑖 𝑋 𝑖 )

= 𝐴0𝑋 0 .𝐵 0𝑋 0 + (𝐴1𝑋 1 .𝐵 0𝑋 0 + 𝐴0𝑋 0 .𝐵 1𝑋 1 )

𝑖=0

+ · · · + (𝐴𝑛−1𝑋 𝑛−1 .𝐵𝑛−2𝑋 𝑛−2 such that 𝑟𝑖 = (𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 ) (mod 𝑋 ), 𝑐 0 = 0 and 𝑐𝑖+1 = ⌊(𝐴𝑖 + 𝐵𝑖 + 𝑐𝑖 )/𝑋 ⌋, and 𝑐 out = 𝑐𝑛+1 □

+ 𝐴𝑛−2𝑋 𝑛−2 .𝐵𝑛−1𝑋 𝑛−1 )

= Proof of correctness of vertical and crosswise multiplication. The DoT-Mul-Words algorithm provides the pseudocode for implementing vertical and crosswise multiplication on SIMD machines. As the registers can accommodate 32- or 64-bit numbers, the results are stored separately in lo and hi registers corresponding to the 128-bit product. We show the correctness of the vertical and crosswise multiplication below:

𝐴=

𝑖

𝐴𝑖 𝑋 = 𝐴𝑛−1𝑋

𝑛−1

+ · · · + 𝐴1 𝑋 + 𝐴0 𝑋

𝑛−1 ∑︁

𝐴𝑖 𝑋 𝑖 𝐵 𝑗 𝑋 𝑗

Thus, DoT-Mul-Words(𝐴, 𝐵, 𝑛, 𝑘) = 𝐴 × 𝐵

□

Probability of carry generation in large-number addition. Lemma B.4 (Probability of a maxed-out limb sum). Let 𝑋𝑖 and 𝑌𝑖 be independent random variables drawn uniformly from {0, 1, . . . , 2𝑘 −1}. The probability that their sum equals the maximum representable value without overflow is Pr 𝑋𝑖 + 𝑌𝑖 = 2𝑘 − 1 = 2−𝑘 . Proof. The sample space is {0, . . . , 2𝑘 −1}2 with 22𝑘 equally likely outcomes. For the event 𝑋𝑖 + 𝑌𝑖 = 2𝑘 −1, fixing any 𝑋𝑖 = 𝑗 uniquely determines 𝑌𝑖 = 2𝑘 −1−𝑗. Since 0 ≤ 𝑗 ≤ 2𝑘 −1 implies 0 ≤ 2𝑘 −1−𝑗 ≤ 2𝑘 −1, this value of 𝑌𝑖 is always valid. There are therefore exactly 2𝑘 favorable pairs, one per value of 𝑋𝑖 , giving

𝑖=0

𝐵=

𝑛−1 ∑︁ 𝑛−1 ∑︁ 𝑗=0 𝑖=0

Theorem B.3 (Correctness of vertical and crosswise multiplication). Suppose two large integers 𝐴 and 𝐵, 𝑛−1 ∑︁

(2)

+ (𝐴𝑛−1𝑋 𝑛−1 .𝐵𝑛−1𝑋 𝑛−1 )

𝐵 𝑗 𝑋 𝑗 = 𝐵𝑛−1𝑋 𝑛−1 + · · · + 𝐵 1𝑋 1 + 𝐵 0𝑋 0

2𝑘 Pr 𝑋𝑖 + 𝑌𝑖 = 2𝑘 − 1 = 2𝑘 = 2−𝑘 . 2

𝑗=0

□

For 𝑘 = 64 this is ≈ 5.4 × 10−20 : a maxed-out limb sum is astronomically rare. Carries themselves, however, are common:

where 𝑋 = 2𝑘 (where 𝑘 is the bit size of a limb), the number of limbs for the two integers is 𝑛, and 0 ≤ 𝐴𝑖 , 𝐵 𝑗 < 𝑋 . Then 𝐴 × 𝐵 = DoT-Mul-Words(𝐴, 𝐵, 𝑛, 𝑘)

Lemma B.5 (Probability of a carry out of a single limb). Under the same conditions,

Proof. Without loss of generality, we assume that both 𝐴 and 𝐵 have the same number of limbs (this can be done by padding the smaller number with required number of 0s on 17

1 Pr 𝑋𝑖 + 𝑌𝑖 ≥ 2𝑘 = − 2− (𝑘+1) . 2

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

8 76

92 81

24 10

□

Bit Size

As 𝑘 → ∞ this approaches 12 : roughly every other limb addition produces a carry. The crucial distinction from Proposition B.4 is that while carries are frequent, carry propagation chains require the intermediate sum 𝑅𝑖 to be maxed out: which is exponentially rare.

Figure 7. Execution time (normalized, lower is better) of DoT (AVX512), two-level KSA (add512/sub512), and Ren et al.’s ProposedAdd/ProposedSub for addition and subtraction across 512– 32768-bit pathological operands.

Corollary B.6 (Phase 4 is negligible for random inputs). Phase 4 of DoT-Add (Algorithm 1) is triggered at limb 𝑖 only when the Phase 1 output 𝑅𝑖 = 2𝑘 −1 and the propagated carry into that limb is 1, so that Phase 3’s addition overflows. By Proposition B.4, Pr[𝑅𝑖 = 2𝑘 −1] = 2−𝑘 . Since limbs are independent, the occurrence of a carry cascade for 𝑛-limb addition is at most (2−𝑘 )𝑛 , which is exponentially small.

2.5x

SSE, 128-bit AVX2, 256-bit

AVX512, 512-bit Baseline (sub-with-bo ow, 64-bit)

0.5x

Figures 7 and 8 show DoT’s performance on pathological inputs, complementing the random-input results in the main text. DoT performs comparably to the two-level KSA and Ren et al.’s method on pathological inputs, confirming that its speedup on random inputs does not come at the cost of worse performance on worst-case inputs. The SIMD speedup plot confirms that AVX512 still achieves a significant speedup over scalar even on pathological inputs, while SSE and AVX2 underperform due to carry-management overhead dominating at narrower widths. DoT’s Impact on Higher-Level Applications. Since OpenSSL’s speed benchmark reports throughput in terms of operations per second, we further compare the latency distributions of DoTSSL and OpenSSL for RSA sign/verify, FFDH derive, and DSA sign/verify (Figure 9). DoTSSL consistently achieves lower latency across all operations. For instance, we observe median (P50) latency improvements reaching up to 7.9% for DSA (2048-bit verify), 6.1% for FFDH derive (2048-bit), and 5.5% for RSA sign (4096-bit). These latency distributions confirm that the throughput gains observed earlier stem directly from faster individual operations, showing consistent performance even at tail percentiles (e.g., 7.8% and 7.7% improvements at the 95th and 99th percentiles for DSA 2048-bit verify).

1.0x

Additional Plots and Results on Intel Xeon Gold 6548Y+ (Emerald Rapids)

1.5x

Speedup

2.0x

2𝑘 − 1 1 Pr 𝑋𝑖 + 𝑌𝑖 ≥ 2𝑘 = 𝑘+1 = − 2− (𝑘+1) . 2 2

𝑘

The carry-producing pairs number 22𝑘 − 2 (22 +1) = 2 (22 −1) . Dividing by 22𝑘 :

Sub [KSA] Add [Ren et a .] Sub [Ren et a .]

2𝑘 (2𝑘 + 1) . 2

𝑠=0

(𝑠 + 1) =

100

Time (ns, Log Sca e)

𝑘 −1 2∑︁

Add [DoT] Sub [DoT] Add [KSA]

Proof. Count non-carry pairs, i.e., those with 𝑋𝑖 + 𝑌𝑖 ≤ 2𝑘 −1. For each value 𝑠 ∈ {0, . . . , 2𝑘 −1}, there are 𝑠 + 1 pairs that sum to exactly 𝑠, so the total number of non-carry pairs is

Bit Size

Figure 8. Speedup of DoT SIMD variants (𝑤=2 SSE, 𝑤=4 AVX2, 𝑤=8 AVX512) over scalar _subborrowx_u64 for subtraction across 512–32768-bit pathological operands.

Results on Intel Xeon Max 9462 (Sapphire Rapids)

This appendix presents the full evaluation results on the Intel Xeon Max 9462 CPU (Sapphire Rapids microarchitecture, SPR), complementing the 6548Y+ (Emerald Rapids, ER) results reported in the main paper. The experimental setup, workloads, and methodology are identical to those described in Section 4. All SPR experiments were compiled with clang 14.0 using -O2 -march=native. D.1

DoT Addition and Subtraction vs. Prior Works

Figure 10(a) shows the execution time comparison of DoT (AVX512) against the two-level KSA and Ren et al.’s method 18

Leveraging SIMD for Accelerating Large-number Arithmetic

1024

2048

3072

4096

Implementation

7680

100% RSA 75% Sign 50% 25% 0% 800

RSA Verify

M 4.5

14M 16M

5.5

Cycles

5M 30M 32.

dot_mul_5 × 5 dot_mul_4 × 4 OpenSSL BN_mul Gueron & Krasnov [37] GMP mpz_mul

M M 180 190

100% 75% 50% 25% 0% K K 110 115

K K K 245 255 380

400

Cycles

M 1.8

1.9

M .3M 5.1 5

M 9.5

10M

Cycles

S-1024

S-2048

V-1024

5M 4.5M 41M 23. 2

221 265 655 350 872

36.0 47.1 86.0 103.5 111.4

6.1 5.6 7.6 3.3 7.8

43M

D.4

V-2048

5M .05M 550K 575K 2

1.9

Cycles

DoT’s Performance over GMP and OpenSSL

Addition and Subtraction. Figure 10(c) shows the speedup of DoTMP over GMP and DoTSSL over OpenSSL on SPR. For addition, DoTMP achieves a 2.80× geomean speedup over GMP (1.84× for smaller operands, 4.25× for larger, peak 4.98×), with instruction count reductions averaging 64.7% (up to 78.5%). DoTSSL achieves 2.92× geomean over OpenSSL (1.92× small, 4.45× large, peak 4.93×), with instruction count reductions of 49.8% on average. For subtraction, DoTMP achieves 3.07× geomean (2.06× small, 4.55× large, peak 5.17×) and DoTSSL achieves 4.00× geomean (2.25× small, 7.10× large, peak 8.24×). These results are lower than ER for DoTMP (ER: 3.81×/3.73×) but comparable for DoTSSL, consistent with SPR’s architectural differences in SIMD throughput.

OpenSSL DoTSSL K K 650 675

5M 5M 1.8 1.9

Figure 9. Latency CDFs of DoTSSL vs. OpenSSL for RSA

sign/verify, FFDH derive, and DSA sign/verify across the evaluated key sizes. Cycles are measured via RDTSC.

on SPR for random test cases. The trends closely mirror those on ER. Compared to the two-level KSA, DoT (AVX512) achieves a geomean speedup of 1.4× for addition (1.23× for smaller operands, 1.73× for larger) and 1.4× for subtraction (1.12× for smaller, 1.73× for larger). Compared to Ren et al., DoT achieves 1.82× geomean speedup for random addition and 1.45× for random subtraction. Instruction count reductions over Ren et al. are 29.1% for random addition and 17.2% for random subtraction, consistent with ER results.

Multiplication. Figure 10(d) shows multiplication speedups on SPR. DoTMP achieves 1.40× geomean over GMP (1.30×– 1.48×, instruction count reduction 47.2% on average). DoTSSL achieves 1.20× geomean over OpenSSL (1.10×–1.32×, instruction count reduction 41.6%). Both are consistent with ER results (1.41× and 1.20×), confirming that multiplication gains are architecture-independent.

DoT’s Performance over SIMD Widths

Figure 10(b) shows the timing speedup of DoT SIMD variants over scalar add-with-carry on SPR. With SSE (𝑤 = 2), geomean speedups are 0.6× for addition and 0.7× for subtraction. With AVX2 (𝑤 = 4), speedups are 1.2× for both. AVX512 (𝑤 = 8) achieves 1.78× for addition and 1.77× for subtraction, slightly below ER’s 1.85×/1.84×, consistent with SPR’s higher base clock and different memory subsystem characteristics. D.3

IPC

bit multiplication on the Intel Xeon Max 9462 (SPR). Cycle counts were obtained via RDTSC/RDTSCP (minimum 700 million total cycles per reported value, averaged over millions of iterations.)

3M 5M 1.4 1.4

100% FFDH 75% Derive 50% 25% 0%

D.2

Avg. Cycles

Table 5. Instruction counts, average cycle counts, and IPC for 256-

40K 50K

100% 75% DSA 50% 25% 0%

Instructions

D.5

DoT’s Impact on Higher-Level Applications

GMPbench. Figure 13 shows DoTMP’s improvement on GMPbench on SPR. The overall score improves by 6.2%. The multiply aggregate improves by 12.7% (individual cases up to 40.0%), divide by 6.9% (up to 30.3%), and pi by 10.1% (up to 13.7% for the 1M-digit case). GCD and GCDext improve by 2.6% and 2.2% respectively. RSA improves by 2.8%. These gains are consistently lower than ER (ER: 7.8% overall) but follow the same workload-dependent pattern.

256-bit Multiplication with DoT

Compared to the ER results, SPR shows a similar pattern of DoT multiplication outperforming both GMP and OpenSSL baselines: Table 4 shows that DoT’s 4 × 4 multiplication routine achieves similar speedups over these baselines on SPR as on ER, confirming that the multiplication gains are architecture-independent and primarily driven by reduced dependency chains rather than memory subsystem differences.

OpenSSL Speed. Figure 14 shows DoTSSL’s throughput improvements on OpenSSL speed on SPR. Improvements are generally higher than on ER: RSA averages 3.9% across sign, verify, encrypt, and decrypt (up to 6.0% for 4096-bit encrypt). FFDH averages 5.4% (up to 7.2% for 4096-bit groups). DSA sign averages 4.4% and verify 5.4% (up to 6.9% for 2048-bit verify). The stronger gains on SPR reflect its higher base 19

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

Method DoT KSA Ren et al. Operation Add Sub

Speedup

Mul (DoTSSL/OpenSSL)

Bit Size

(d) Mul in integrated libs

6 57

96 40

2 51

8 76

6 57

57 24

38 16

2 51

8 76

1.2x

38 16

Mul (DoTMP/GMP)

1.4x

1.0x

0x 15

Bit Size

(b) Add SIMD variants speedup

Operation Add Sub

1.0x

Bit Size

1.5x

(a) Add/sub vs prior works

2.0x

0.5x

Pair DoTMP/GMP DoTSSL/OpenSSL

SSE (w=2) AVX2 (w=4) AVX512 (w=8)

92 81

96 40

Speedup

100

Time (ns) (log scale)

2.5x

Figure 10. Micro-benchmark evaluation of DoT on the Intel Xeon Max 9462 (SPR). (a) Execution time (log scale) of DoT (AVX512), two-level

Add [DoT] Sub [DoT] Add [KSA]

100

Sub [KSA] Add [Ren et a .] Sub [Ren et a .]

2.5x

SSE, 128-bit AVX2, 256-bit

AVX512, 512-bit Baseline (sub-with-bo ow, 64-bit)

2.0x

Speedup

Time (ns, Log Sca e)

KSA, and Ren et al. across 512–32768-bit random operands. (b) Speedup of DoT SIMD variants (𝑤=2 SSE, 𝑤=4 AVX2, 𝑤=8 AVX512) over scalar _addcarryx_u64 for addition. (c) Timing speedup of DoTMP over GMP and DoTSSL over OpenSSL for addition and subtraction. (d) Timing speedup of DoTMP over GMP and DoTSSL over OpenSSL for multiplication.

1.5x 1.0x

Bit Size

8 76

6 57

4 38

8 28

2 51

8 76 32

4 38 16

92 81

96 40

48 20

24 10

0.5x

Bit Size

Figure 11. Execution time (normalized, lower is better) of DoT (AVX512), two-level KSA, and Ren et al.’s method for addition and subtraction across 512–32768-bit pathological operands on the Intel Xeon Max 9462 (SPR), pathological test cases.

Figure 12. Speedup of DoT SIMD variants (𝑤=2 SSE, 𝑤=4 AVX2, 𝑤=8 AVX512) over scalar _subborrowx_u64 for subtraction across 512–32768-bit pathological operands on the Intel Xeon Max 9462 (SPR).

frequency making the relative cost of scalar carries more pronounced. Similarly, the latency distributions (Figure 16) show that DoTSSL consistently achieves lower latency across all operations on SPR. For instance, we observe median (P50) latency

improvements reaching up to 6.1% for DSA (2048-bit verify), 6.7% for FFDH derive (4096-bit), and 5.8% for RSA sign (4096-bit). These latency improvements maintain consistent performance even at tail percentiles (e.g., 6.2% and 6.1% improvements at the 95th and 99th percentiles for DSA 2048-bit verify). 20

Leveraging SIMD for Accelerating Large-number Arithmetic

<1.5% -2.1% +6.5% +12.4% +18.4% -1.7%

128 512 8K 128K 2048K 128x128 512x512 Multiply 8Kx8K 128Kx128K 2048Kx2048K 15Kx10K 20Kx10K 30Kx10K 16384Kx512 16384Kx256K

GCD

128K

128

+40.0% +25.5% +16.5% +16.3% +9.1% +16.8%

512

GCDext

<1.5% <1.5% -2.7%

Pi 100K

−5

(a) Mul/Div Improvement (%)

+6.9% +12.7%

app

+6.4%

base

+6.0%

total

+6.2%

+13.7%

1000K

+12.9% 10

+7.7% +9.1%

10K

+30.3%

+2.6%

Multiply

+9.9%

+16.0%

GCD

+4.2% <1.5% +3.9%

512

<1.5%

+2.8% +2.2%

+4.9%

1024K

RSA

RSA GCDext

Divide

8K 128K

<1.5% <1.5% <1.5% <1.5%

+11.2%

1024K

<1.5% +6.0%

−5

<1.5% +2.9%

+10.1%

512

+37.0%

8K÷32 8K÷64 8K÷128 8K÷4K Divide128K÷64K 8192K÷4096K 8K÷8064 16384K÷256K

<1.5% <1.5%

128

−5

(b) GCD/RSA/Pi Improvement (%)

Figure 13. DoTMP’s percentage improvement over GMP across GMPbench workloads on the Intel Xeon Max 9462 (SPR). Overall score

improves by 6.2%, with multiply (+12.7%) and pi (+10.1%) leading, following the same workload-dependent pattern as ER but at modestly lower absolute gains.

DoT’s Contribution to these gains. Similar to ER, we used perf to analyze the cycle composition of DoT’s routines in GMPbench and OpenSSL speed on SPR. Figure 15 shows that DoT’s dot_add_words, dot_sub_words, and dot_mul_4x4

Sign/s Verify/s

Improvement (%)

Encrypt/Encaps Decrypt/Decaps

routines account for a significant share of cycles in both GMPbench and OpenSSL speed workloads, confirming that DoT’s improvements in these core operations are driving the overall performance gains observed in higher-level applications on SPR as well. Encrypt/Encaps Decrypt/Decaps

Keygen (op/s)

0 1024

2048

3072

4096

Key Size (bits) (a) RSA

7680

0 1024

2048

3072

4096

Key Size (bits) (b) RSA KEM

7680

Sign/s Verify/s

0 2048

3072

4096

6144

Group Size (bits) (c) FFDH

8192

1024

2048

Key Size (bits) (d) DSA

Figure 14. DoTSSL throughput improvement (%) over OpenSSL for RSA, RSA KEM, FFDH, and DSA on the Intel Xeon Max 9462 (SPR). Improvements are generally higher than on ER: FFDH reaches up to +7.2% and DSA verify up to +6.9%, reflecting SPR’s higher base frequency amplifying the relative cost of scalar carry chains. 21

Subhrajit Das, Abhishek Bichhawat, and Yuvraj Patel

512

GCD

512×512

1024

128K

2048

RSA 3072

4096

8K×8K 15K×10K

GCDext

20K×10K

7680

128K 1M

Multiply30K×10K 128K 128K×128K

512

RSA

2M×2M

DSA

1024 2048

16M×512

dot_mul_4x4 dot_add_words dot_sub_words

16M×256K 128K÷64K

10K

2048

Divide 8M÷4M

Pi 100K

FFDH 3072

16M÷256K

4096

1M 0

Cycles Spent (%) in DoT Routines

(a) GMPbench (Mul, Div)

(b) GMPbench (GCD, RSA, Pi)

Figure 15. Cycle spent (%) by DoT’s dot_add_words, dot_sub_words, and dot_mul_4x4 routines in GMPbench and OpenSSL speed workloads, measured via perf on the Intel Xeon Max 9462 (SPR). We omitted handful of cases in the GMPbench (e.g., lower sized mul, div and gcd) since they spend zero cycles in DoT routines. Additionally, OpenSSL speed benchmarks keygen, sign, encrypt, decrypt, etc. in aggregate for each key size; thus we report the cycle share for DoT routines in the aggregate of all operations for each key size. 1024 RSA Sign

2048

3072

4096

7680

100% 75% 50% 25% 0% 1M .25M 1

18M 20M 35M 40M

Cycles

220

240

100% RSA 75% Verify 50% 25% 0% 40K

50K

K K 140 150

K K 310 320

K K 500 525

1.8

Cycles

FFDH Derive

5 1.8

100% 75% 50% 25% 0% M M 2.3 2.4

100% 75% DSA 50% 25% 0%

S-1024

K K 800 850

M 75M 12M 6.

6.5

S-2048

14M

Cycles

V-1024

30M 32M 52.5M

5M 57.

V-2048 OpenSSL DoTSSL

M M 2.5 2.75

K K 700 725

M M 2.4 2.5

Cycles

Figure 16. Latency comparison (CDF) of DoTSSL vs. OpenSSL for RSA sign/verify, FFDH derive, and DSA sign/verify. Cycles are measured via RDTSC on the Intel Xeon Max 9462 (SPR) and plotted on a log scale.

Record · ID 126482 · SHA-256 deb2086d03ffa060
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.