Vectorized SQIsign Implementation Using AVX-512

Vectorized SQIsign Implementation Using AVX-512 Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

arXiv:2608.13948v1 [cs.CR] 14 Aug 2026

Fudan University, Shanghai, China, {wzwang23,ctwang24,wuyu25,qfxue24,jyzheng23}@m.fudan.edu.cn,ylzhao@fudan.edu.cn

Abstract. SQIsign is the sole isogeny-based digital signature scheme submitted to the NIST Post-Quantum Cryptography standardization process, distinguished by its foundation on the hardness of the endomorphism ring problem for supersingular elliptic curves. While offering compact key and signature sizes, SQIsign’s practical deployment is hindered by computationally intensive signing procedures. This paper presents the first comprehensive vectorized implementation of SQIsign utilizing the AVX-512 Integer Fused Multiply-Add (IFMA) instruction set architecture. By systematically redesigning the computational stack—encompassing prime-field and extension-field arithmetic, elliptic curve operations including batched point doubling and scalar multiplication, as well as pairing computations via cubical arithmetic and two-dimensional isogeny evaluations—we achieve substantial performance improvements over the reference implementation. When combined with Qlapoti technology, our implementation attains a 2.69× speedup for signing and a 3.18× improvement for verification at NIST security level I. Contrary to misconceptions regarding the obsolescence of AVX-512, we emphasize that Intel’s AVX10 instruction set architecture (revision 10.2, scheduled for widespread deployment in late 2026) will standardize AVX-512 capabilities—including IFMA instructions—across both performance and efficiency cores, ensuring long-term viability of these optimization techniques. Furthermore, our vectorization strategies are architecture-agnostic and provide a methodological foundation applicable to broader isogeny-based cryptographic constructions. This work demonstrates that SIMD vectorization represents a critical yet underexplored optimization dimension for postquantum isogeny-based schemes, independent of recent algorithmic advances. Keywords: SQIsign · Vectorization · Isogeny · AVX-512

1

Introduction

The advent of quantum computing poses a fundamental threat to classical public-key cryptography, particularly to schemes based on integer factorization and discrete logarithm problems. In response to this challenge, the National Institute of Standards and Technology (NIST) initiated a standardization process for post-quantum cryptographic algorithms, focusing on schemes that resist attacks by both classical and quantum computers. Among the candidates for digital signatures, SQIsign [AAA+ 24] stands out as the sole isogeny-based submission, distinguished by its foundation on the hardness of the endomorphism ring problem for supersingular elliptic curves. While offering compact key and signature sizes, SQIsign instantiates a Fiat-Shamir transform [FS87] over a computationally intensive sigma protocol that translates quaternion algebra ideals into explicit isogeny computations [BDD+ 24], necessitating extensive arithmetic over large prime fields

2

Vectorized SQIsign Implementation Using AVX-512

Fp and their quadratic extensions Fp2 , elliptic curve operations, and higher-dimensional isogeny evaluations over principally polarized abelian surfaces (PPAS). SQIsign instantiates a Fiat-Shamir transform [FS87] over a sigma protocol proving knowledge of an isogeny path between supersingular curves. The signature generation process involves three computationally intensive phases: commitment, challenge, and response. The performance bottleneck resides primarily in the IdealToIsogeny conversion routine, which translates quaternion algebra ideals into explicit isogeny computations [BDD+ 24]. This routine necessitates extensive arithmetic over large prime fields Fp and their quadratic extensions Fp2 , elliptic curve operations including point doubling and scalar multiplication, and higher-dimensional isogeny evaluations over principally polarized abelian surfaces (PPAS). Consequently, optimizing these low-level arithmetic primitives is essential for rendering SQIsign practical for real-world deployment. To address these computational demands, we leverage the Advanced Vector Extensions 512 (AVX-512) instruction set architecture, specifically the Integer Fused MultiplyAdd (IFMA) extension. AVX-512 represents Intel’s fourth-generation SIMD architecture, expanding the x86-64 execution environment to encompass thirty-two 512-bit vector registers, effectively doubling both the register count and vector width compared to AVX2. The IFMA extension, available since the Ice Lake microarchitecture and broadly supported in server-grade Xeon Scalable processors as well as AMD’s Zen 5 microarchitecture, provides specialized instructions vpmadd52luq and vpmadd52huq that execute eight parallel 52-bit × 52-bit multiplications with fused multiply-add capabilities. These instructions are particularly adept at accelerating multi-precision modular arithmetic, eliminating separate carry propagation instructions and significantly reducing instruction count during Montgomery multiplication. Importantly, AVX-512 is not a deprecated technology: Intel has announced the AVX10 instruction set architecture, with revision 10.2 scheduled for widespread deployment in late 2026, which will standardize AVX-512 capabilities— including IFMA instructions—across both performance (P) cores and efficiency (E) cores, ensuring universal availability of these cryptographic acceleration primitives in future processor generations [Int26]. Contributions. This work presents the first comprehensive vectorized implementation of SQIsign utilizing the AVX-512IFMA instruction set, achieving substantial performance improvements over the reference implementation. Our contributions are threefold: 1. Comprehensive Vectorization of SQIsign. We have systematically redesigned the computationally intensive components of SQIsign to exploit vectorized finite field arithmetic across all NIST security levels (I, III, and V). By vectorizing the majority of time-consuming operations in the elliptic curve, pairing, and isogeny computation layers, our implementation achieves speedups ranging from 1.03× to 1.28× for key generation, 1.08× to 1.26× for signing, and 1.46× to 1.84× for verification compared to the reference implementation. 2. Extensive Vectorization Analysis Beyond ARM Implementations. De Feo et al. [FJWY26] demonstrated the feasibility of SQIsign vectorization on ARM architectures using NEON instructions, achieving significant speedups through vectorization of the IdealToIsogeny subroutine. Our work explores complementary vectorization strategies on the x86-64 platform with AVX-512, encompassing a significantly broader range of functions including elliptic curve scalar multiplication using Montgomery ladders, batched point doubling on Jacobian coordinates, bilinear pairing computations via cubical arithmetic, and theta-coordinate operations for two-dimensional isogenies. These vectorization techniques are architecture-agnostic in their algorithmic design, allowing their migration to other SIMD platforms including future ARM SVE or RISC-V vector extensions.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

3

3. Generalizable Optimizations for Isogeny-Based Cryptography. Our optimizations for elliptic curve operations—including parallel xDBLADD implementations, batched scalar multiplication strategies, and vectorized theta-point arithmetic— are not restricted to SQIsign. The techniques we develop for point doubling in various coordinate systems, simultaneous evaluation of pairings, and isogeny computations provide a methodological foundation applicable to other isogeny-based cryptosystems utilizing higher-dimensional abelian varieties. Related Works. The optimization of isogeny-based cryptography through vectorization has been explored in several contexts. Cheng et al. [CFGR22] demonstrated highly vectorized implementations of SIKE using AVX-512IFMA, establishing the foundation for our finite field arithmetic vectorization as well as isogeny operation vectorization. Their work on Montgomery multiplication and reduction serves as a basis for our Fp and Fp2 arithmetic implementations. De Feo et al. [FJWY26] presented the first vectorized implementation of SQIsign on ARM Cortex-A76 architectures, achieving significant speedups through NEON-based vectorization of the IdealToIsogeny subroutine. Our work independently explores vectorization opportunities on the x86-64 platform with AVX-512, providing analysis of additional computational components including elliptic curve pairings and theta-coordinate operations. Furthermore, the recent work of Kim et al. [KLKL25] on fixed-precision integer arithmetic for SQIsign eliminates dependencies on the GMP library. While the present work focuses on the “curve-side” computations (finite fields, elliptic curves, and isogenies), the combination of vectorized field arithmetic with fixed-precision quaternion operations represents a promising avenue for future research and comprehensive end-to-end optimization. Outline. The remainder of this paper is organized as follows. Section 2 reviews the mathematical foundations of SQIsign, including the Deuring correspondence and the protocol workflow, alongside a detailed description of the AVX-512IFMA architecture and its relevance to cryptographic implementations. Section 3 details our vectorized implementations of prime-field Fp and quadratic extension-field Fp2 arithmetic. Section 4 presents vectorized algorithms for elliptic curve operations, including point doubling in various coordinate systems, scalar multiplication, and pairing operations based on cubical arithmetic. Section 5 addresses the vectorization of isogeny computations in both dimension one and dimension two, covering the algorithmic adaptations required for efficient SIMD execution. Section 6 presents our experimental results and performance analysis, comparing our AVX-512 implementation against the reference implementation. Finally, Section 7 concludes the paper and outlines future research directions.

2

Preliminaries

2.1

SQIsign

SQIsign represents a promising candidate in the ongoing NIST post-quantum digital signature standardization process, distinguished by its reliance on the hardness of the endomorphism ring problem for supersingular elliptic curves—a mathematical foundation fundamentally distinct from lattice-based and hash-based alternatives. At its core, the scheme instantiates a Fiat-Shamir transform over a sigma protocol proving knowledge of an isogeny path between supersingular curves. Mathematical Foundation. Let p be a prime of cryptographic size, and consider the finite field Fp2 . A supersingular elliptic curve E defined over Fp2 possesses the property

4

Vectorized SQIsign Implementation Using AVX-512

that its endomorphism ring End(E) is isomorphic to a maximal order in a quaternion algebra Bp,∞ ramified at p and infinity. The Deuring correspondence establishes an equivalence between isogenies originating from a fixed base curve E0 and left ideals of the corresponding maximal order O0 ∼ = End(E0 ). This correspondence enables the translation of algebraic operations in quaternion algebras into geometric operations on elliptic curves, and vice versa. Protocol Overview. The signature generation process comprises three phases: commitment, challenge, and response. The prover generates a secret isogeny φsk : E0 → Epk whose codomain constitutes the public key. During signing, the prover samples a random commitment isogeny φcom : E0 → Ecom . Upon receiving a challenge isogeny φchl : Epk → Echl from the verifier, the prover must demonstrate knowledge of a response isogeny φrsp : Ecom → Echl that does not backtrack through φchl . Crucially, the prover utilizes the Deuring correspondence to compute an ideal representation of the composite isogeny φchl ◦ φsk ◦ φ̂com , then samples an equivalent ideal amenable to efficient isogeny translation. Computational Components. The performance characteristics of SQIsign are dominated by the Ideal-to-Isogeny conversion routine, which comprises two distinct computational components. The first involves arithmetic over quaternion algebras—specifically, ideal equivalence testing, lattice reduction, and norm computation—which is currently implemented using the GNU Multiple Precision Arithmetic Library (GMP) for arbitraryprecision integer arithmetic. Recent work by Kim et al. [KLKL25] demonstrates that these quaternion operations can be reformulated using fixed-precision integer arithmetic, eliminating the dependency on GMP and laying the groundwork for future SIMD vectorization. However, we do not optimize this algebraic component in the present work, leaving it as a potential target for subsequent vectorization efforts. The second component—and the primary focus of this work—comprises the “curveside” computations: finite field arithmetic in Fp and Fp2 , elliptic curve operations (point doubling, scalar multiplication), and higher-dimensional isogeny evaluations over principally polarized abelian surfaces (PPAS). Contemporary instantiations of SQIsign employ two-dimensional isogeny chains derived from Kani’s lemma [Kan97], necessitating extensive computations on Jacobian varieties and theta coordinates. The translation between ideal representations and these explicit geometric objects represents the primary bottleneck amenable to vectorization, as these operations ultimately reduce to modular multiplication and reduction in large prime fields. Implementation Scope. Figure 1 illustrates the modular dependency structure of the SQIsign implementation. We employ the following color coding to indicate our optimization scope: red components comprise the quaternion algebra operations relying on GMP, along with the fixed-precision integer module present in the current SQIsign implementation, neither of which are optimized in this work; green components denote algorithms that have been redesigned for vectorized execution; and blue components represent higher-level operations that derive performance benefits from their dependencies on our vectorized field arithmetic primitives and vectorized algorithms.

2.2

AVX-512

The Advanced Vector Extensions 512 (AVX-512) constitute Intel’s fourth-generation SIMD architecture, introducing substantial enhancements over preceding AVX and AVX2 generations. These extensions provide the computational infrastructure necessary for high-

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

5

Figure 1: SQIsigns dependency graph from [FJWY26]. throughput, low-latency implementations of cryptographic primitives operating on large integer operands. Architectural Overview. AVX-512 expands the x86-64 execution environment to encompass thirty-two 512-bit vector registers (ZMM0-ZMM31), effectively doubling both the register count and vector width compared to AVX2. The architecture supports operations on packed data elements ranging from 8-bit integers to 64-bit integers, enabling simultaneous processing of either eight 64-bit or sixteen 32-bit quantities within a single instruction. The AVX-512 Foundation (AVX-512F) provides the core integer and floatingpoint capabilities, while various sub-extensions address specific computational domains. Integer Fused Multiply-Add (IFMA). Of particular relevance to cryptographic implementations is the AVX-512IFMA extension, introduced with the Cannon Lake microarchitecture and broadly available in Ice Lake and subsequent server-grade processors. This extension furnishes two specialized instructions: vpmadd52luq (Packed Multiply-Add 52-bit Unsigned, Lower half) and vpmadd52huq (Packed Multiply-Add 52-bit Unsigned, Higher half).These instructions execute eight parallel 52-bit × 52-bit multiplications, producing eight 104-bit intermediate products. The vpmadd52luq instruction adds the lower 52 bits of each product to the corresponding 64-bit destination element, while vpmadd52huq adds the upper 52 bits. This fused multiply-add capability eliminates separate addition instructions for carry propagation in multi-precision arithmetic, significantly reducing instruction count and port pressure during Montgomery multiplication and modular reduction operations. Hardware Support and Future Prospects. While AVX-512IFMA support remains limited in consumer desktop processorsprimarily appearing in 10th generation (Ice Lake) and 11th generation (Tiger Lake) Core seriesthe extension enjoys broad support in servergrade Xeon Scalable processors and AMD’s Zen 5 microarchitecture, which implements full 512-bit vector datapaths [Adv24]. Looking ahead, Intel has announced the AVX10 instruction set architecture, with revision 10.2 scheduled for widespread deployment in late 2026. This converged vector ISA will standardize AVX-512 capabilities, including IFMA instructions, across both performance (P) cores and efficiency (E) cores, ensuring universal availability of these cryptographic acceleration primitives in future processor generations[Int26].

3

Finite Field Arithmetic

The performance of isogeny-based cryptographic schemes, including SQIsign, is fundamentally constrained by arithmetic in the underlying finite fields. This section describes

6

Vectorized SQIsign Implementation Using AVX-512

our AVX-512 implementation of the prime-field Fp and quadratic extension-field Fp2 operations. Our approach follows the vectorization strategies established in prior work on isogeny-based cryptography [CFGR22], utilizing the (8 × 1)-way parallelization pattern at the Fp layer.

3.1

Representation and Data Layout

The arithmetic framework of our implementation builds upon Montgomery’s modular multiplication methodology [Mon85], which eliminates expensive division operations through a precomputed radix-dependent inverse. We adopt a radix-251 representation for field elements, motivated by the interplay between the 52-bit multiplier width of AVX-512IFMA and the requirement for lazy reduction. Field elements are decomposed into five 51-bit limbs, stored within the lower 52 bits of each 64-bit lane to prevent overflow during IFMA operations. The use of 51-bit limbs (rather than the full 52 bits) reserves headroom for carry accumulation during consecutive operations, thereby enabling lazy reduction strategies at higher algorithmic layers. 8-Way Limb Vector Organization. Our vectorized Fp implementation adopts an (8 × 1)way parallelization strategy, processing eight independent field elements simultaneously within a single 512-bit ZMM register. Following the approach in [CFGR22], we assign each 64-bit lane to a distinct field element, maintaining spatial correspondence between the logical and physical representation. Taking the NIST Level I parameter set as an example, given eight field elements a(0) , . . . , a(7) ∈ Fp , each represented as five 51-bit limbs, the limb vector set V is organized as:  (0) (1) (2) (3) (4) (5) (6) (7)  [a0 , a0 , a0 , a0 , a0 , a0 , a0 , a0 ]     (0) (1) (2) (3) (4) (5) (6) (7)      [a1 , a1 , a1 , a1 , a1 , a1 , a1 , a1 ] (1) (2) (3) (4) (5) (6) (7) V = [a(0) , a , a , a , a , a , a , a ] 2 2 2 2 2 2 2 2    (0) (1) (2) (3) (4) (5) (6) (7)   [a3 , a3 , a3 , a3 , a3 , a3 , a3 , a3 ]      (0) (1) (2) (3) (4) (5) (6) (7)  [a4 , a4 , a4 , a4 , a4 , a4 , a4 , a4 ] where each row constitutes a 512-bit ZMM register containing the i-th limb of all eight elements.

3.2

Prime-Field Fp Operations

The special form of SQIsign’s characteristic prime p = c · 2f − 1 (with small c) renders it amenable to optimized Montgomery reduction where the multiplicative factor involves only the small constant c [FLOR18]. This Montgomery-friendly property simplifies the modular reduction step. Multiplication and Reduction. The modular multiplication leverages the special form p = c · 2f − 1 to streamline the reduction phase. For SQIsign’s parameters, the Montgomery reduction of a 510-bit intermediate product proceeds via multiplication by the precomputed constant, implemented using a sequence of vpmadd52luq and vpmadd52huq instructions. The five-limb structure of our operands ensures that the reduction step requires only a minimal number of vector operations, with carries propagated across the 51-bit boundaries through mask-shift-add sequences.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

3.3

7

Quadratic Extension Fp2 Layer

Operations in Fp2 are built upon the Fp primitives using the quadratic extension with irreducible polynomial X 2 + 1 (exploiting the fact that p ≡ 3 (mod 4)). We implement three distinct parallelization strategies for Fp2 multiplication and squaring, corresponding to (8×1×1)-way, (4×2×1)-way, and (2×4×1)-way execution patterns, where the notation (w × x × y) indicates w parallel Fp2 operations, each utilizing x parallel Fp operations, with y denoting the lane utilization at the lowest level. • (8 × 1 × 1)-way: Eight independent Fp2 multiplications proceed in parallel, with each Fp operation internally utilizing the (8 × 1)-way limb organization. This mode maximizes throughput for batch operations. • (4×2×1)-way: Four Fp2 multiplications execute concurrently, with each multiplication internally processing two Fp operations (e.g., real and imaginary components) in parallel. This hybrid approach balances parallelism with register pressure for intermediate computations. • (2 × 4 × 1)-way: Two Fp2 multiplications proceed with four-fold internal parallelism, suitable for algorithms requiring simultaneous processing of multiple field elements within a single higher-level operation. The choice among these parallelization modes depends on the specific requirements of the higher-level algorithm. The (8 × 1 × 1)-way approach excels when processing independent batches of field elements, such as in simultaneous evaluation of isogenies at multiple points. The (4 × 2 × 1)-way and (2 × 4 × 1)-way variants are preferred when the algorithmic structure demands internal parallelism within individual Fp2 operations, such as in point addition formulas where real and imaginary components can be computed concurrently. Each implementation employs either schoolbook or Karatsuba decomposition at the coefficient level, selected based on the trade-off between multiplication count and addition/subtraction overhead under the specific parallelization constraints.

3.4

Addition and Subtraction

Vectorized implementation of Fp2 addition and subtraction is straightforward. The addition r = a + b = (a0 + b0 ) + (a1 + b1 )i in Fp2 consists of two independent additions in Fp , which are executed in parallel using the vpaddq instruction. Similarly, subtraction employs the vpsubq instruction. AVX-512 provides mask registers (k-registers) that enable selective operation on specific lanes of a ZMM register. By utilizing masked vpaddq and vpsubq instructions, we can perform addition on certain lanes while performing subtraction on others, deferring the modular correction until a unified reduction step. This capability allows us to mix addition and subtraction operations within the same algorithmic step (e.g., computing both a + b and a − b simultaneously, or handling mixed sequences of additions and subtractions that arise in algorithms such as the Hadamard operator) without requiring intermediate reductions or sign corrections.

3.5

Performance Evaluation

Table 1 presents the cycle counts for Fp2 multiplication and squaring across different implementation variants, measured under the NIST Level I parameter set (p = 5 · 2248 − 1). The comparison includes: (1) hand-optimized x64 assembly utilizing 64-bit limbs with mulx/adcx/adox instructions; (2) C-based implementation generated by the modarith framework using 51-bit limbs but without vectorization; and (3) our AVX-512IFMA vectorized implementations with varying degrees of parallelism at the Fp layer.

8

Vectorized SQIsign Implementation Using AVX-512

Table 1: Performance comparison of Fp2 arithmetic implementations for NIST Level I (cycles/operation). Operation Mul Sqr

x64 ASM 102 72

modarith 246 158

AVX-512IFMA 8-way

4-way

2-way

38 22

46 25

51 49

The results demonstrate that the vectorized implementation achieves substantial speedups over both scalar variants. Notably, the performance gap between the 51-bit non-vectorized C implementation and the 64-bit assembly baseline stems primarily from the difference in limb count: for the 248-bit prime p, the 64-bit representation requires only 4 limbs, whereas the 51-bit representation requires 5 limbs. This difference in limb count directly impacts the number of elementary multiplications required in the finite field multiplication algorithms.

3.6

Imperative for Comprehensive Vectorization

Modern x86 microarchitectures provide powerful scalar instructions for multi-precision arithmetic, specifically the mulx, adcx, and adox instructions available since the Broadwell generation. These instructions enable highly efficient 64-bit limb implementations with optimized carry propagation, rendering non-vectorized 51-bit limb approaches—even with algorithmic optimizations—at a significant performance disadvantage. Consequently, to realize the full potential of the 51-bit representation within the AVX512IFMA framework, it is insufficient to vectorize only isolated field operations. While not all operations in SQIsign can be efficiently vectorized due to data dependencies and other factors, the majority of computations on elliptic curves and PPAS have been restructured to operate on vectorized data. Frequent transitions between scalar 51-bit and vectorized representations incur substantial overhead from data format conversion and register spilling, necessitating a careful redesign of the algorithmic stack to minimize such transitions and maintain data in vector registers throughout the bulk of the computation.

4

Elliptic Curve Arithmetic

In this section, we review the fundamental elliptic curve operations underlying SQIsign, including point doubling, scalar multiplication, pairings, and discrete logarithm computations. For efficiency in isogeny arithmetic, SQIsign predominantly employs Montgomery curves of the form BY 2 Z = X 3 + AX 2 Z + XZ 2 equipped with projective x-only coordinates (X : Z). The curve parameter A is also represented projectively as A = (APr : CPr ) with A = APr /CPr , and A24 is represented projectively as (A24 : C24 ) = (APr + 2CPr : 4CPr ). However, certain high-dimensional isogeny computations require access to the full (x, y) coordinates of points. In such cases, SQIsign transitions to Jacobian coordinates on Montgomery curves, or alternatively, to Modified Jacobian coordinates[CMO98] on short Weierstrass curves y 2 = x3 +ax+b. The coordinate transformations are defined as follows: • For Jacobian coordinates (X : Y : Z): x = X/Z 2 , y = Y /Z 3 • For Modified Jacobian coordinates (X : Y : Z : T ): x = X/Z 2 , y = Y /Z 3 , T = aZ 4 The remainder of this section is organized as follows. Subsection 4.1 presents vectorized implementations of point doubling across various coordinate models. Subsection 4.2 turns

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

9

to scalar multiplication, describing vectorized algorithms for both the Montgomery ladder and two-dimensional multi-scalar multiplication. Subsection 4.3 addresses the vectorization of pairing operations employed in elliptic curve discrete logarithm computations.

4.1

Point Doubling

We begin with point doubling in x-only coordinates. Given the curve parameter (A24 : C24 ) , a single doubling operation requires 2 squarings and 4 multiplications. Due to data dependencies, at most two multiplications or squarings can be executed simultaneously. In the reference implementation of SQIsign, a normalization procedure is invoked when more than 50 consecutive doublings are performed, reducing C24 to 1 and thereby decreasing the multiplication count per doubling to 3. However, this optimization offers no benefit for our vectorized implementation; consequently, we consistently omit the normalization step when performing vectorized point doublings. Algorithm 1 Vectorized xDBL_A24 Input: A projective point P = (XP : ZP ) and the Montgomery coefficient (A24 : C24 ) of the curve E. Output: The projective point [2]P = (X2P : Z2P ). s0 ← XP − ZP 1: t0 ← XP + ZP s1 ← s20 2: t1 ← t20 3: t2 ← t1 − s1 s3 ← B × C24 4: t3 ← t2 × A24 5: t4 ← t3 + s3 s5 ← t2 × t4 6: t5 ← s3 × t1 7: return [2]P = (t5 : s5 ) During the signing and verification procedures of SQIsign, multiple steps require simultaneously doubling a full-order basis to adjust them to an appropriate order. This entails performing point doublings on three x-only coordinate points on the same curve, an operation evidently more amenable to vectorization due to its higher degree of parallelism. Without normalization, doubling three points requires 6 squarings and 12 multiplications in total; after normalization, which reduces C24 to 1, the multiplication count decreases to 9. This reduction translates to a decrease from 5 to 4 multiplication rounds under fourway parallelization. Given that the number of doublings required for order adjustment of bases is typically substantial, we consistently perform the normalization step prior to executing the operation. The overwhelming majority of elliptic curve operations in SQIsign are performed using x-only coordinates on Montgomery curves. However, during higher-dimensional isogeny computations, complete projective coordinate arithmetic becomes occasionally necessary. The most common scenario arises in the computation of (2, 2)-isogeny chains, where iterative doublings of two points are required (see [AAA+ 24, Algorithm 8.46]). In the current SQIsign implementation, such operations are mostly accomplished by mapping points to modified Jacobian coordinates on short Weierstrass curves, with each point doubling incurring a cost of 3 finite field multiplications and 5 finite field squarings. In Algorithm 3, we present a simultaneous doubling of two points with the multiplications and squarings organized into maximally balanced vectorized operations. In [FJWY26], De Feo et al. analyze this setting and argue that the aforementioned formula is ill-suited for vectorization. They instead propose an alternative formula requiring 4 finite field multiplications and 4 finite field squarings, upon which they design their vectorization strategy. This approach appears to demand fewer field additions and subtractions compared to our method. We defer a comprehensive performance comparison

10

Vectorized SQIsign Implementation Using AVX-512

Algorithm 2 Doubling three points on the same curve. Input: Projective points P = (XP : ZP ), Q = (XQ : ZQ ), R = (XR : ZR ) and the Montgomery coefficient (A24 : C24 ) of the curve E. Output: The projective point [2]P = (X2P : Z2P ), [2]Q = (X2Q : Z2Q ), [2]R = (X2R : Z2R ). s0 ← XP − ZP m0 ← XQ + ZQ n0 ← XQ − ZQ 1: t0 ← XP + ZP s1 ← s20 m1 ← m20 n1 ← n20 2: t1 ← t20 s2 ← m1 − n1 m2 ← XR + ZR n2 ← XR − ZR 3: t2 ← t1 − s1 s3 ← s2 × A24 m3 ← m2 × m2 n3 ← n2 × n2 4: t3 ← t2 × A24 s4 ← s3 + n1 m4 ← m3 − n3 5: t4 ← t3 + s1 s5 ← s4 × s2 m5 ← m4 × A24 6: t5 ← t4 × t2 7: t6 ← m5 + n3 s7 ← m1 × n1 m7 ← m3 × n3 n7 ← m4 × t6 8: t7 ← t1 × s1 9: return [2]P = (t7 : t5 ), [2]Q = (s7 : s5 ), [2]R = (m7 : n7 )

to subsequent work. Algorithm 3 Vectorized batched DBLW Input: Modified Jacobian points P = (XP : YP : ZP : TP ), Q = (XQ : YQ : ZQ : TQ ) for a Weierstrass curve. Output: The doubled modified Jacobian points [2]P, [2]Q. 2 2 s0 ← XQ m0 ← YP2 n0 ← YP2 1: t0 ← XP s1 ← s0 + s0 m1 ← m0 + m0 n1 ← n0 + n0 2: t1 ← t0 + t0 s2 ← s1 + s0 m2 ← m1 + XP n2 ← n1 + XQ 3: t2 ← t1 + t0 s3 ← n22 m3 ← m21 n3 ← n21 4: t3 ← m22 s4 ← s2 + TQ m4 ← m3 + m3 n4 ← n3 + n3 5: t4 ← t2 + TP s5 ← s3 − s0 6: t5 ← t3 − t0 s6 ← s5 − n3 7: t6 ← t5 − m3 s7 ← s4 × s4 m7 ← YP × ZP n7 ← YQ × ZQ 8: t7 ← t4 × t4 s8 ← s7 − s6 9: t8 ← t7 − t6 s9 ← s8 − n6 10: t9 ← t8 − m6 s10 ← s6 − n9 11: t10 ← t6 − t9 s11 ← s4 × s10 m11 ← m4 × TP n11 ← n4 × TQ 12: t11 ← t4 × t10 s12 ← s11 − n4 13: t12 ← t11 − m4 s13 ← n7 + n7 m13 ← m11 + m11 n13 ← n11 + n11 14: t13 ← m7 + m7 15: return [2]P = (t9 : t12 : t13 : m13 ), [2]Q = (s9 : s12 : s13 : n13 )

4.2

Scalar Multiplication

Scalar multiplication constitutes a fundamental operation employed across various subroutines within the SQIsign signature scheme. The SQIsign specification categorizes scalar multiplication algorithms into three distinct classes, all of which operate on Montgomery curves utilizing projective x-only coordinate representations: Ladder The classical Montgomery ladder algorithm, originally introduced by Montgomery in [Mon87]. Given a point P on an elliptic curve E, this algorithm computes the scalar multiple [m]P . Ladder3pt The three-point ladder algorithm proposed by Faz-Hernández et al. in [FLOR18]. Given points P , Q, and P − Q on an elliptic curve E, this algorithm computes P + [m]Q.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

11

LadderBiscalar The two-dimensional ladder algorithm due to Bernstein [Ber06]. Given points P , Q, and P − Q on an elliptic curve E, this algorithm computes the biscalar multiple [m]P + [n]Q. For the sake of brevity, we refrain from presenting the detailed theoretical foundations underlying these algorithms. Instead, we excerpt the (partial) pseudocode specifications as documented in the SQIsign reference implementation. Algorithm 4 Ladder(P, E, m) Input: A projective point P = (XP : ZP ), the Montgomery coefficient (A24 : C24 ) of the curve E, and a positive scalar m with binary representation m = (mk1 , . . . , m0 )2 . Output: The projective point [m]P = (X[m]P : Z[m]P ). 1: ((X0 : Z0 ), (X1 : Z1 )) ← ((1 : 0), (XP : ZP )) 2: for i from k − 1 down to 0 do 3: if mi = 1 then 4: ((X1 , Z1 ), (X0 , Z0 )) ← xDBLADD((X1 : Z1 ), (X0 : Z0 ), (XP : ZP ), (A24 : C24 )) 5: else 6: ((X0 : Z0 ), (X1 : Z1 )) ← xDBLADD((X0 : Z0 ), (X1 : Z1 ), (XP : ZP ), (A24 : C24 )) 7: (X[m]P : Z[m]P ) ← (X0 : Z0 ) 8: return [m]P = (X[m]P : Z[m]P )

Algorithm 5 Ladder3pt(P, Q, P − Q, (A24 : C24 ), m) Input: Projective points P = (XP : ZP ), Q = (XQ : ZQ ), P − Q = (XP −Q : ZP −Q ), the Montgomery coefficient (A24 : C24 ) of the curve E, and a positive scalar m with binary representation m = (mk1 , . . . , m0 )2 . Output: The projective point P + [m]Q = (XP +[m]Q : ZP +[m]Q ). 1: ((X0 : Z0 ), (X1 : Z1 ), (X2 : Z2 )) ← ((XP : ZP ), (XQ : ZQ ), (XP −Q : ZP −Q )) 2: for i from 0 up to k − 1 do 3: if mi = 1 then 4: ((X0 , Z1 ), (X0 , Z0 )) ← xDBLADD((X0 : Z0 ), (X1 : Z1 ), (X2 : Z2 ), (A24 : C24 )) 5: else 6: ((X0 : Z0 ), (X2 : Z2 )) ← xDBLADD((X0 : Z0 ), (X2 : Z2 ), (X1 : Z1 ), (A24 : C24 )) 7: (XP +[m]Q : ZP +[m]Q ) ← (X1 : Z1 ) 8: return P + [m]Q = (XP +[m]Q : ZP +[m]Q ) A careful examination of the Ladder and Ladder3pt algorithms reveals that their main loops exhibit strict data iteration dependencies, preventing parallel unrolling at the loop level; consequently, their computational performance is entirely determined by the execution efficiency of the underlying xDBLADD function.. Therefore, the key to enhancing algorithmic efficiency through vectorization lies in exploiting the internal parallelism within the xDBLADD operation. The xDBLADD function constitutes a combined operation that simultaneously performs point doubling and point addition. It takes as input the curve parameters together with points P , Q, and P −Q, and produces as output the values P +Q and [2]P . In the literature [CFGR22], Cheng et al. conducted a comparative analysis of four-way vectorization strategies for the xDBLADD procedure as described in [CS09, NS22, HEY22], concluding that the vectorization approaches presented in [CS09] and [NS22] are heavily dependent on the special form of the Montgomery curve parameter A24 —a characteristic commonly encountered in classical elliptic curve cryptography—rendering them unsuitable for scalar multiplication implementations in isogeny-based cryptography. By contrast, the strategy proposed in [HEY22] was deemed more appropriate for this context.

12

Vectorized SQIsign Implementation Using AVX-512

XQ

ZQ

XP

H

×

H

×

×

×

−

H

×

ZP

×

× A24 × C24 +

× ZP −Q × XP −Q XP +Q

ZP +Q

×

×

X2P

Z2P

Figure 2: Our vectorized strategy for xDBLADD However, the algorithm in [HEY22] relies on the curve parameter A, whereas in isogenybased cryptography, curve parameters are typically expressed in the form (A : C) or (A24 : C24 ) for computational efficiency. To address this limitation, we present a redesigned four-way vectorization strategy for xDBLADD in Algorithm 6. This implementation imposes no special requirements on the parameter format and eliminates the need for parameter normalization operations; the parameters need only be provided in the (A24 : C24 ) representation. To facilitate comparison with [NS22] and [HEY22], we additionally provide a graphical illustration of the algorithm in Figure 2. Algorithm 6 xDBLADD algorithm with 4-way parallelization at Fp2 -level. Input: Projective points P = (XP : ZP ), Q = (XQ : ZQ ), P − Q = (XP −Q : ZP −Q ), and the Montgomery coefficient (A24 : C24 ) of the curve E. Output: Projective points [2]P = (X2P : Z2P ) and P + Q = (XP +Q : ZP +Q ). s0 ← XQ − ZQ m0 ← XP + ZP n0 ← XP − ZP 1: t0 ← XQ + ZQ s1 ← s0 × m0 m1 ← m0 × m0 n1 ← n0 × n0 2: t1 ← t0 × n0 s2 ← t1 − s1 m2 ← m1 − n1 3: t2 ← t1 + s1 s3 ← s2 × s2 m3 ← m2 × A24 n3 ← n1 × C24 4: t3 ← t2 × t2 5: t4 ← m3 + n3 s5 ← s3 × XP −Q m5 ← m1 × n3 n5 ← m2 × t4 6: t5 ← t3 × ZP −Q
7: return ([2]P, P + Q) = (m5 : n5 ), (t5 : s5 ) We now proceed to analyze the LadderBiscalar algorithm. Owing to its substantially greater complexity compared to the aforementioned ladder algorithms, we present in Algorithm 7 only the core portion subject to our primary optimizations; the complete algorithmic specification may be found in Algorithm 8.8 in [AAA+ 24]. To facilitate comprehension, we begin with a brief overview of the algorithmic work-

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

13

flow. The fundamental distinction between this algorithm and its predecessors lies in the scalar processing mechanism: whereas previous ladder algorithms process input scalars bit-by-bit, with each loop iteration governed solely by the current bit value, LadderBiscalar exhibits significantly more intricate behavior. The operations performed in each loop iteration depend not merely on the “current bit”, but additionally on the operation executed in the preceding iteration as well as the subsequent bit. This renders LadderBiscalar a “stateful” algorithm in a meaningful sense. Consequently, the algorithm naturally decomposes into two distinct phases: the Recoding stage and the Evaluation stage. During the Recoding stage, the algorithm scans the input scalars m and n to determine the specific operation to be executed in each loop iteration. Subsequently, the Evaluation stage commences with the elliptic curve computations. The initialization procedure employs the points (1 : 0), P , and Q—in a specific order—to establish the values R0 , R1 , and R2 , while simultaneously initializing four difference points D1 , D2 , F1 , and F2 . Throughout the iterative process, the values R0 , R1 , and R2 undergo transformation, whereas the difference points remain invariant. The computational workload of each loop iteration consists of selecting one point from {R0 , R1 , R2 } for point doubling, together with two differential additions performed on pairs of points. This computation encompasses a total of 18 finite field multiplications. In Algorithm 8, we reorganize this main loop into five four-way parallel multiplication operations. As normalization of parameters offers no reduction in the number of parallel multiplications, we maintain our design choice of parameter non-normalization. Algorithm 7 LadderBiscalar(P, Q, P − Q, (A24 : C24 ), m, n) (Partial) Input: Projective points P = (XP : ZP ), Q = (XQ : ZQ ), P − Q = (XP −Q : ZP −Q ), the Montgomery coefficient (A24 : C24 ) of the curve E, and positive scalars m, n with binary representations m = (mk−1 , . . . , m0 )2 and n = (nk−1 , . . . , n0 )2 , resp. Output: The projective point [m]P + [n]Q = (X[m]P +[n]Q : Z[m]P +[n]Q ). 29: for i from k − 1 down to 0 do 30: h ← r2i + r2i+1 31: T0 ← Rh (mod 2) 32: T = (T0 , T1 ) ← (T0 , R2 ) 33: T0 ← xDBL(T⌊h/2⌋ , (A24 : C24 )) 34: T1 ← Rr2i+1 35: T2 ← Rr2i+1 +1 36: if r2i+1 = 1 then 37: T M P ← D1 38: D1 ← D2 39: D2 ← T M P 40: T1 ← xADD(T1 , T2 , D1 ) 41: T2 ← xADD(R0 , R2 , F1 ) 42: if h (mod 2) = 1 then 43: T M P ← F1 44: F1 ← F2 45: F2 ← T M P 46: R0 ← T0 47: R1 ← T1 48: R2 ← T2 49: (X[m]P +[n]Q : Z[m]P +[n]Q ) ← R(m (mod 2)⊕1)+(n (mod 2)⊕1) 50: return [m]P + [n]Q = (X[m]P +[n]Q : Z[m]P +[n]Q )

14

Vectorized SQIsign Implementation Using AVX-512

Algorithm 8 LadderBiscalar main loop algorithm with 4-way parallelization at Fp2 -level. Input: Projective points (XT0 : ZT0 ), (XT1 : ZT1 ), (XT2 : ZT2 ), (XR0 : ZR0 ), (XR2 : ZR2 ), (XD1 : ZD1 ), (XF1 : ZF1 ). Output: Projective points R0 = [2]T0 , R1 = T1 + T2 , R2 = R0 + R2 . s0 ← XT2 − ZT2 m0 ← XT1 + ZT1 n0 ← XT1 − ZT1 1: t0 ← XT2 + ZT2 s1 ← XR0 − ZR0 m1 ← XR2 + ZR2 n1 ← XR2 − ZR2 2: t1 ← XR0 + ZR0 s2 ← t0 × n0 m2 ← t1 × n1 n2 ← s1 × m1 3: t2 ← s0 × m0 s3 ← t2 − s2 m3 ← m2 + n2 n3 ← m2 − n2 4: t3 ← t2 + s2 s4 ← s3 × s3 m4 ← m3 × m3 n4 ← n3 × n3 5: t4 ← t3 × t3 s5 ← XT0 − ZT0 6: t5 ← XT0 + ZT0 s6 ← s4 × XD1 m6 ← t5 × t5 n6 ← s5 × s5 7: t6 ← t4 × ZD1 8: t7 ← m6 − n6 s8 ← n6 × C24 m8 ← m4 × ZF1 n8 ← n4 × XF1 9: t8 ← t7 × A24 10: t9 ← t8 + s8 s10 ← t7 × t9 11: t10 ← m6 × s8 12: return R0 = (t10 : s10 ), R1 = (t6 : s6 ), R2 = (m8 : n8 )

4.3

Pairing and Discrete Logarithm

To accelerate computations, the key generation and signing procedures in SQIsign require the generation of a 2e -torsion basis (Q1 , Q2 ) on a specific elliptic curve, which is subsequently recorded within the private key or signature. To conserve bandwidth, SQIsign does not store the coordinates of the basis points directly. Instead, the scheme proceeds from a deterministic basis (P1 , P2 ) of E[2e ] and computes the change-of-basis matrix      x1 x2 P1 R = . x3 x4 P2 S This matrix could theoretically be generated via discrete logarithm methods on the elliptic curve; however, for computational efficiency, SQIsign adopts an alternative approach involving the computation of five pairings followed by discrete logarithms in the finite field. We excerpt this algorithm in Algorithm 9. Algorithm 9 ChangeOfBasis2e (E, (P1 , P2 ), (Q1 , Q2 )) Input: A basis (P1 , P2 ) and a basis (Q1 , Q2 ) for E[2e ]. Output: A change-of-basis matrix (xi ),with 1 ≤ i ≤ 4, so that Q1 = [x1 ]P1 + [x2 ]P2 and Q2 = [x3 ]P1 + [x4 ]P2 . 1: ζ ← t2e (P1 , P2 ) 2: ζ1 ← t2e (Q1 , P2 ), ζ2 ← 1/t2e (Q1 , P1 ), ζ3 ← t2e (Q2 , P2 ), ζ4 ← 1/t2e (Q2 , P1 ) 3: for i from 1 up to 4 do 4: xi ← 2f −e · logζ (ζi ) 5: return (x1 , x2 , x3 , x4 )

The current implementation of SQIsign employs the Tate-Lichtenbaum pairing to facilitate discrete logarithm computations. Notably, the pairing computation in SQIsign does not follow the classical Miller loop algorithm; instead, it utilizes an approach based on cubical arithmetic. The theoretical foundations of cubical arithmetic trace back to reinterpretations of the Tate pairing as monodromy in biextensions associated to principal polarizations, as developed in [Sta08] and subsequently refined in [Rob24] through the lens of cubical arithmetic structures.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

15

We refrain from presenting the complete theoretical framework of cubical arithmetic herein, confining our discussion to the core algorithmic steps relevant to SQIsign. The implementation employs level-2 cubical arithmetic associated to the divisor 2(0E ). To e above P and Q. This inicompute tn (P, Q), one first obtains cubical points Pe and Q tialization is straightforward: the points are normalized to yield Pe = (X(P )/Z(P ) : 1) e = (X(Q)/Z(Q) : 1). Subsequently, the CubicalLadder algorithm is executed to and Q e The CubicalLadder algorithm constitutes essencompute the cubical points 2e Pe and 2e Q. tially a cubical arithmetic adaptation of the three-point ladder, wherein CubicalDiffAdd and CubicalDbl denote the cubical arithmetic variants of differential addition and point doubling, respectively. The detailed specifications appear in [AAA+ 24, Algorithms 8.13 and 8.14]. Algorithm 10 CubicalLadder(E, e, P^ + Q, Pe, x(Q)) Input: An elliptic curve E : y 2 = x3 + Ax2 + x in Montgomery model, an integer e, two cubical points P^ + Q = (X(P + Q), Z(P + Q)), Pe = (X(P ), Z(P )) and the xcoordinate x(Q) of Q. e Output: The cubical points 2e Pe and 2e Pe + Q. ^ 1: nP Q ← P +Q e 2: nP ← P 3: for k from 1 up to e do 4: nP Q ← CubicalDiffAdd(E, nP Q, nP, x(Q)) 5: nP ← CubicalDbl(E, nP ) 6: return (nP, nP Q) Both differential addition and point doubling in cubical arithmetic require five finite field multiplications, which we reorganize into three four-way parallel vectorized multiplication operations. Pope et al. [PRR+ 25] observed that SQIsign invariably computes two bilinear pairings with respect to the same point (cf. Algorithm 9, Line 2), thereby enabling the reuse of intermediate values. Capitalizing on this insight, we fuse one point doubling with two differential additions to obtain Algorithm 11, which requires three fourway parallel multiplications and one four-way parallel squaring operation. The resulting ladder implementation can simultaneously generate the two bilinear pairings required in Algorithm 9. Algorithm 11 cubicalDBLADDADD e = (XQ , ZQ ), R e = (XR , ZR ), the inverses Input: Three cubical points Pe = (XP , ZP ), Q of the differences ixP R = 1/x(P − R), ixQR = 1/x(Q − R) and the normalized Montgomery coefficient A24 of the curve E. e and the cubical differential additions P^ ^ Output: The cubical double 2R + R,Q +R 1: t1 ← XP + ZP s1 ← XP − ZP m1 ← XR + ZR n1 ← XR − ZR 2: t2 ← XQ + ZQ s2 ← XQ − ZQ 3: t3 ← m1 × m1 s3 ← n1 × n1 m3 ← m1 × s2 n3 ← n1 × t2 4: t4 ← t3 − s3 s4 ← m3 + n3 m4 ← m3 − n3 5: t5 ← t1 × n1 s5 ← s1 × m1 m5 ← A24 × t4 6: t6 ← t5 + s5 s6 ← t5 − s5 m6 ← s3 + m5 7: t7 ← t6 × t6 s7 ← s6 × s6 m7 ← s4 × s4 n7 ← m4 × m4 8: t8 ← ixP R × t7 s8 ← t3 × s3 m8 ← m6 × t4 n8 ← ixQR × m7 e ^ ^ 9: return 2R = (s8 : m8 ),P + R = (t8 : s7 ),Q + R = (n8 : n7 )

16

Vectorized SQIsign Implementation Using AVX-512

Remark. The execution of the CubicalDbl algorithm necessitates normalized curve parameters A24 . One might observe that CubicalDbl is in fact identical to the x-only coordinate doubling formula on Montgomery curves, and consequently surmise that the normalization operation could be eliminated by applying our previously established technique— introducing a multiplication by C24 within Algorithm 11—without increasing the number of vectorized multiplications while saving one field inversion. This approach, however, proves untenable. Such a modification to the doubling formula alters the normalization of e 0, which in turn necessitates corresponding adjustments to the cubical differential addition formula; continued use of the original formulae would yield erroneous results. A comprehensive analysis of this phenomenon is provided in [Rob24, Section 5.2].

5

Isogenies in Dimension 1 and 2

5.1

Isogenies between Elliptic Curves

The evaluation of 4-isogenies on Montgomery curves follows the formulas established by Costello et al. [CLN16] and subsequently optimized for vectorization by Cheng et al. [CFGR22]. Given a kernel point P4 = (X4 : Z4 ) of order 4 and a point P = (XP : ZP ) on the domain curve, the image point P ′ = ϕ4 (P ) = (XP ′ : ZP ′ ) is computed via the relations:   XP ′ = 16 · (X4 XP − Z4 ZP )2 + Z42 (XP2 − ZP2 ) · (X4 XP − Z4 ZP )2 ,   ZP ′ = 16 · (X4 ZP − Z4 XP )2 − Z42 (XP2 − ZP2 ) · (X4 ZP − Z4 XP )2 . Algebraically, these formulas expose at most two independent multiplications that can execute simultaneously—specifically, the computation of X4 XP ±Z4 ZP and X4 ZP ±Z4 XP followed by their squares. This structural constraint limits the theoretical parallelism per point evaluation to two-way, as recognized in prior work [CFGR22]. Cheng et al. [CFGR22] observed that isogeny chain computations using optimal strategy traversal naturally maintain a stack of points awaiting evaluation at each step. Rather than processing these points sequentially with repeated scalar-to-vector data marshalling between steps, they demonstrated that batching multiple evaluations into a single vectorized invocation amortizes instruction overhead and improves throughput. We adopt this batching methodology in our SQIsign implementation, interleaving coordinates from multiple points across AVX-512 lanes to saturate the SIMD multipliers. Table 3 quantifies this effect, demonstrating that processing multiple points concurrently yields speedups, despite the two-way structural constraint inherent to the formulas.

5.2

Isogenies between Principally Polarized Abelian Surfaces

SQIsign’s two-dimensional isogeny computations operate on theta coordinates, requiring intensive point doubling operations on the PPAS. The theta model, originating from Mumford’s seminal work on algebraic theta functions [Mum83], provides efficient arithmetic on principally polarized abelian varieties through level-2 theta constants. Lubicz and Robert [LR12] subsequently developed computational frameworks for isogeny evaluation within this coordinate system, which we adapt for our vectorized implementation. We provide vectorized implementations for both the precomputation phase (Algorithm 12) and the doubling operation (Algorithm 13), reorganizing the Hadamard transforms and coordinate multiplications to maximize SIMD lane utilization. For the generic evaluation of (2, 2)-isogenies, our derived formulas align with the vectorized approach presented in the ARM implementation [FJWY26, Algorithm 16]. Following the methodology established in the one-dimensional case, we apply the same batching principle to two-dimensional isogeny evaluations, processing multiple theta points concurrently

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

Algorithm 12 ThetaPrecomp_vec(0A ) Input: Theta null point 0A := (a : b : c : d). Output: Auxiliary constants consts used for arithmetic. s0 ← b2 m0 ← c2 1: t0 ← a2 s1 ← t0 − s0 m1 ← m0 + n0 2: t1 ← t0 + s0 s2 ← t1 − m1 m2 ← s1 + n1 3: t2 ← t1 + m1 s3 ← c · d m3 ← t2 · m2 4: t3 ← s2 · n2 abd ← n3 · d acd ← s3 · a 5: abc ← n3 · c ABD ← m3 · n2 ACD ← t3 · t2 6: ABC ← m3 · s2 7: consts← {abc, abd, acd, bcd, ABC, ABD, ACD, BCD} 8: return consts

17

n0 ← d2 n1 ← m0 − n0 n2 ← s1 − n1 n3 ← a · b bcd ← s3 · b BCD ← t3 · m2

Algorithm 13 ThetaDBL_vec(P ,consts) Input: The theta coordinates of P on A with theta null point 0A := (a : b : c : d), and the auxiliary constants consts:=ThetaPrecomp(0A ). Output: The theta coordinates of the point [2]P . 1: xP , yP , zP , wP ← P 2: c1 , c2 , c3 , c4 , c5 , c6 , c7 , c8 ← consts s0 ← yP2 m0 ← zP2 n0 ← wP2 3: t0 ← x2P s1 ← t0 − s0 m1 ← m0 + n0 n1 ← m0 − n0 4: t1 ← t0 + s0 s2 ← t1 − m1 m2 ← s1 + n1 n2 ← s1 − n1 5: t2 ← t1 + m1 s3 ← m22 m3 ← s22 n3 ← n22 6: t3 ← t22 s4 ← s3 · c7 m4 ← m3 · c6 n 4 ← n 3 · c5 7: t4 ← t3 · c8 s5 ← t4 − s4 m5 ← m4 + n4 n5 ← m4 − n4 8: t5 ← t4 + s4 s6 ← t5 − m5 m6 ← s5 + n5 n6 ← s5 − n5 9: t6 ← t5 + m5 Y2P ← m6 · c3 Z2P ← s6 · c2 W2P ← n6 · c1 10: X2P ← t6 · c4 11: return (X2P : Y2P : Z2P : W2P )

to amortize data movement costs. The performance characteristics for these operations are summarized in Table 4.

6

Experimental Results

We implemented our vectorized SQIsign using AVX-512IFMA atop the official reference implementation1 . All benchmarks were conducted on an Intel Core i7-11700F CPU (2.50GHz base frequency) running Ubuntu 24.04 LTS with GCC 13.2.0, using the project’s default release configuration with hyper-threading and Turbo Boost disabled. Our AVX-512IFMA implementation uniformly employs a radix-251 representation for finite field elements. Since the modarith2 generator produces 51-bit limbs for NIST Level I but divergent limb sizes for Levels III and V, we modified the scripts to generate 51-bit limb implementations across all parameter sets, ensuring compatible data layouts between vectorized and scalar code paths. For comparison, the “ASM” entries in the tables denote the Broadwell-optimized reference implementation from the official SQIsign library, which employs 64-bit limbs via hand-written assembly for finite field arithmetic. Table 2 presents the performance of core elliptic curve operations across NIST security levels I, III, and V. The “Iter.” column indicates the iteration count specific to each operation type: doublings per point for xDBL, and scalar bit length for ladder operations. 1 https://github.com/SQIsign/the-sqisign 2 https://github.com/mcarrickscott/modarith

18

Vectorized SQIsign Implementation Using AVX-512

Notably, these iteration values are selected based on typical operational requirements for each parameter set rather than fixed constants, reflecting the actual computational loads encountered during signature generation and verification. Tables 3 and 4 demonstrate the throughput scaling of one-dimensional and two-dimensional isogeny evaluations under varying degrees of vectorization. As expected, increased parallelism generally correlates with improved performance, validating our batched evaluation strategy. Table 5 reports the end-to-end performance of our AVX-512IFMA implementation, including the variant combined with the Qlapoti[BCE+ 25] technique. We observe speedups across all three operations—key generation, signing, and verification—for all NIST parameter sets. The verification procedure exhibits particularly substantial improvements. This disparity arises from the differing computational compositions: key generation and signing incorporate quaternion algebra operations that are currently executed via GMP without vectorization, whereas verification relies exclusively on elliptic curve arithmetic and PPAS operations. Consequently, verification benefits maximally from our vectorized finite field and isogeny primitives, free from the serialization bottlenecks imposed by multi-precision integer libraries.

7

Conclusions

This paper presented the first comprehensive vectorization of SQIsign utilizing AVX512IFMA, demonstrating that SIMD optimization represents a critical yet previously underexploited dimension for isogeny-based signature schemes. By systematically redesigning the finite field arithmetic layer and extending vectorization through elliptic curve operations, pairings, and both one- and two-dimensional isogenies, we achieved substantial performance gains across all NIST security levels. The orthogonal nature of our optimizations to recent algorithmic advances suggests that future SQIsign implementations may compound these gains through synergistic application of multiple improvement dimensions. Specifically, our vectorization techniques operate independently of high-level algorithmic refinements such as Qlapoti [BCE+ 25], Superglue [Dup25]. This independence implies that the speedups reported in this work represent a multiplicative factor atop gains achievable through protocol-level innovations. Furthermore, our methodology demonstrates architectural portability that transcends specific SQIsign variants. Whether applied to the original KLPT-based construction [DKL+ 20], the high-dimensional approach of SQIsignHD [DLRW24], the variants SQIsign2DEast [NOC+ 24] or SQIPrime [DF24], the fundamental vectorization strategies for finite field and isogeny arithmetic remain equally applicable. This universality stems from the fact that all these constructions ultimately rely on the same underlying primitives— modular multiplication over Fp2 , theta-coordinate operations, and elliptic curve arithmetic— that constitute our optimization targets. Several avenues remain for future work. The integration of fixed-precision quaternion arithmetic following Kim et al. [KLKL25] could enable vectorization of the IdealToIsogeny algebraic components currently executing in scalar GMP, potentially eliminating the performance gap between signing and verification observed in our experiments. Additionally, the adaptation of our vectorization techniques to ARM SVE or emerging RISC-V vector extensions would validate the architectural portability of our approach. As isogeny-based cryptography continues to mature, we expect vectorization to play an increasingly central role in bridging the gap between theoretical compactness and practical deployment efficiency.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

19

Acknowledgements The authors would like to thank the anonymous reviewers of TCHES (IACR Transactions on Cryptographic Hardware and Embedded Systems), Volume 2026, Issue 4, for their careful reading of this manuscript and their valuable comments and constructive suggestions, which have helped improve this work.

References [AAA+ 24] Marius A. Aardal, Gora Adj, Diego F. Aranha, Andrea Basso, Isaac Andrés Canales Martínez, Jorge Chávez-Saab, Maria Corte-Real Santos, Pierrick Dartois, Luca De Feo, Max Duparc, Jonathan Komada Eriksen, Tako Boris Fouotsa, Décio Luiz Gazzoni Filho, Basil Hess, David Kohel, Antonin Leroux, Patrick Longa, Luciano Maino, Michael Meyer, Kohei Nakagawa, Hiroshi Onuki, Lorenz Panny, Sikhar Patranabis, Christophe Petit, Giacomo Pope, Krijn Reijnders, Damien Robert, Francisco Rodríguez Henríquez, Sina Schaeffler, and Benjamin Wesolowski. SQIsign. Technical report, National Institute of Standards and Technology, 2024. available at https://csrc.nist. gov/Projects/pqc-dig-sig/round-2-additional-signatures. [Adv24]

Advanced Micro Devices. Software Optimization Guide for the AMD Zen5 Microarchitecture, revision 1.0 edition, August 2024. Document Number: 58455.

[BCE+ 25] Giacomo Borin, Maria Corte-Real Santos, Jonathan Komada Eriksen, Riccardo Invernizzi, Marzio Mula, Sina Schaeffler, and Frederik Vercauteren. sfQlapoti: Simple and efficient translation of quaternion ideals to isogenies. In Goichiro Hanaoka and Bo-Yin Yang, editors, ASIACRYPT 2025, Part IV, volume 16248 of LNCS, pages 174–205. Springer, Singapore, December 2025. [BDD+ 24] Andrea Basso, Luca De Feo, Pierrick Dartois, Antonin Leroux, Luciano Maino, Giacomo Pope, Damien Robert, and Benjamin Wesolowski. SQIsign2D-West: The fast, the small, and the safer. Cryptology ePrint Archive, Report 2024/760, 2024. [Ber06]

Daniel J. Bernstein. Differential addition chains. [Online], 2006. https: //cr.yp.to/ecdh/diffchain-20060219.pdf.

[CFGR22] Hao Cheng, Georgios Fotiadis, Johann Großschädl, and Peter Y. A. Ryan. Highly vectorized SIKE for AVX-512. IACR TCHES, 2022(2):41–68, 2022. [CLN16]

Craig Costello, Patrick Longa, and Michael Naehrig. Efficient algorithms for supersingular isogeny Diffie-Hellman. In Matthew Robshaw and Jonathan Katz, editors, CRYPTO 2016, Part I, volume 9814 of LNCS, pages 572–601. Springer, Berlin, Heidelberg, August 2016.

[CMO98]

Henri Cohen, Atsuko Miyaji, and Takatoshi Ono. Efficient elliptic curve exponentiation using mixed coordinates. In Kazuo Ohta and Dingyi Pei, editors, ASIACRYPT’98, volume 1514 of LNCS, pages 51–65. Springer, Berlin, Heidelberg, October 1998.

[CS09]

Neil Costigan and Peter Schwabe. Fast elliptic-curve cryptography on the cell broadband engine. In Bart Preneel, editor, AFRICACRYPT 09, volume 5580 of LNCS, pages 368–385. Springer, Berlin, Heidelberg, June 2009.

20

[DF24]

Vectorized SQIsign Implementation Using AVX-512

Max Duparc and Tako Boris Fouotsa. SQIPrime: A dimension 2 variant of SQISignHD with non-smooth challenge isogenies. In Kai-Min Chung and Yu Sasaki, editors, ASIACRYPT 2024, Part III, volume 15486 of LNCS, pages 396–429. Springer, Singapore, December 2024.

[DKL+ 20] Luca De Feo, David Kohel, Antonin Leroux, Christophe Petit, and Benjamin Wesolowski. SQISign: Compact post-quantum signatures from quaternions and isogenies. In Shiho Moriai and Huaxiong Wang, editors, ASIACRYPT 2020, Part I, volume 12491 of LNCS, pages 64–93. Springer, Cham, December 2020. [DLRW24] Pierrick Dartois, Antonin Leroux, Damien Robert, and Benjamin Wesolowski. SQIsignHD: New dimensions in cryptography. In Marc Joye and Gregor Leander, editors, EUROCRYPT 2024, Part I, volume 14651 of LNCS, pages 3–32. Springer, Cham, May 2024. [Dup25]

Max Duparc. Superglue: Fast formulae for (2,2)-gluing isogenies. In Goichiro Hanaoka and Bo-Yin Yang, editors, ASIACRYPT 2025, Part IV, volume 16248 of LNCS, pages 372–400. Springer, Singapore, December 2025.

[FJWY26] Luca De Feo, Li-Jie Jian, Ting-Yuan Wang, and Bo-Yin Yang. SQISign on ARM. Cryptology ePrint Archive, Paper 2026/394, 2026. [FLOR18] Armando Faz-Hernández, Julio Cesar López-Hernández, Eduardo OchoaJiménez, and Francisco Rodríguez-Henríquez. A faster software implementation of the supersingular isogeny diffie-hellman key exchange protocol. IEEE Transactions on Computers, 67(11):1622–1636, 2018. [FS87]

Amos Fiat and Adi Shamir. How to prove yourself: Practical solutions to identification and signature problems. In Andrew M. Odlyzko, editor, Advances in Cryptology — CRYPTO’ 86, pages 186–194, Berlin, Heidelberg, 1987. Springer Berlin Heidelberg.

[HEY22]

Hüseyin Hışıl, Berkan Eğrice, and Mert Yassı. Fast 4 way vectorized ladder for the complete set of montgomery curves. International Journal of Information Security Science, 11(2):12–24, 2022.

[Int26]

Intel Corporation. Intel® Advanced Vector Extensions 10.2 Architecture Specification, revision 6.0 edition, January 2026. Document Number: 361050006US.

[Kan97]

Ernst Kani. The number of curves of genus two with elliptic differentials. Journal für die reine und angewandte Mathematik, 1997(485):93–122, 1997.

[KLKL25] Won Kim, Jeonghwan Lee, Hyeonhak Kim, and Changmin Lee. SQIsign with fixed-precision integer arithmetic. Cryptology ePrint Archive, Report 2025/1649, 2025. [LR12]

David Lubicz and Damien Robert. Computing isogenies between abelian varieties. Compositio Mathematica, 148(5):14831515, 2012.

[Mon85]

Peter L. Montgomery. Modular multiplication without trial division. Mathematics of Computation, 44(170):519–521, 1985.

[Mon87]

Peter L. Montgomery. Speeding the pollard and elliptic curve methods of factorization. Mathematics of Computation, 48(177):243–264, 1987.

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

[Mum83]

21

David Mumford. Tata Lectures on Theta II, volume 43 of Progress in Mathematics. Birkhäuser, 1983.

[NOC+ 24] Kohei Nakagawa, Hiroshi Onuki, Wouter Castryck, Mingjie Chen, Riccardo Invernizzi, Gioella Lorenzon, and Frederik Vercauteren. SQIsign2D-East: A new signature scheme using 2-dimensional isogenies. In Kai-Min Chung and Yu Sasaki, editors, ASIACRYPT 2024, Part III, volume 15486 of LNCS, pages 272–303. Springer, Singapore, December 2024. [NS22]

Kaushik Nath and Palash Sarkar. Efficient 4-way vectorizations of the montgomery ladder. IEEE Transactions on Computers, 71(3):712–723, 2022.

[PRR+ 25] Giacomo Pope, Krijn Reijnders, Damien Robert, Alessandro Sferlazza, and Benjamin Smith. Simpler and faster pairings from the montgomery ladder. Cryptology ePrint Archive, Report 2025/672, 2025. [Rob24]

Damien Robert. Fast pairings via biextensions and cubical arithmetic. Cryptology ePrint Archive, Report 2024/517, 2024.

[Sta08]

Katherine E. Stange. Elliptic nets and elliptic curves. PhD thesis, Brown University, May 2008.

22

Vectorized SQIsign Implementation Using AVX-512

Table 2: Elliptic curve operation performance (point doublings, ladders) by NIST level. Algorithm

ec_dbl_iter

ec_dbl_iter_basis

DBLW

ec_mul

ec_ladder3pt

ec_biscalar_mul

Level

Iter.

I

120

III

182

V

245

I

120

III

182

V

245

I

123

III

186

V

248

I

128

III

194

V

255

I

256

III

384

V

512

I

248

III

376

V

500

Impl.

Cycles

Speedup

ASM AVX-512 ASM AVX-512 ASM AVX-512

65198 59046 222310 157548 425704 289502

1.00× 1.10× 1.00× 1.41× 1.00× 1.47×

ASM AVX-512 ASM AVX-512 ASM AVX-512

195574 108278 666676 317378 1279354 602438

1.00× 1.81× 1.00× 2.10× 1.00× 2.12×

ASM AVX-512 ASM AVX-512 ASM AVX-512

240406 117354 748552 328336 1409952 602044

1.00× 2.05× 1.00× 2.28× 1.00× 2.34×

ASM AVX-512 ASM AVX-512 ASM AVX-512

146224 94644 526508 277738 982576 512496

1.00× 1.54× 1.00× 1.90× 1.00× 1.92×

ASM AVX-512 ASM AVX-512 ASM AVX-512

297622 193676 1073598 563822 2049516 1049416

1.00× 1.54× 1.00× 1.90× 1.00× 1.95×

ASM AVX-512 ASM AVX-512 ASM AVX-512

499962 321898 1717500 925660 3265876 1743904

1.00× 1.55× 1.00× 1.86× 1.00× 1.87×

Weize Wang, Chutong Wang, Yu Wu, Qifan Xue, Jieyu Zheng and Yunlei Zhao

23

Table 3: 4-isogeny evaluation throughput under varying parallelization degrees. Level

Impl.

#Inst.

Cycles

Cyc./Inst.

Speedup

I

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

828 646 1054 1530 2786

828 646 527 383 348

1.00× 1.28× 1.57× 2.16× 2.38×

III

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

2026 1216 1986 2980 5572

2026 1216 993 745 697

1.00× 1.67× 2.04× 2.72× 2.91×

V

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

2914 1654 2798 4152 7856

2914 1654 1399 1038 982

1.00× 1.76× 2.08× 2.81× 2.97×

Table 4: (2, 2)-isogeny evaluation performance with different vectorization strategies. Level

Impl.

#Inst.

Cycles

Cyc./Inst.

Speedup

I

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

938 654 1176 1894 3232

938 654 588 474 404

1.00× 1.43× 1.60× 1.98× 2.32×

III

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

2072 1128 2326 3694 6244

2072 1128 1163 924 781

1.00× 1.84× 1.78× 2.24× 2.65×

V

ASM AVX-512 AVX-512 AVX-512 AVX-512

1 1 2 4 8

2890 1492 3238 5058 8750

2890 1492 1619 1265 1094

1.00× 1.94× 1.79× 2.29× 2.64×

24

Vectorized SQIsign Implementation Using AVX-512

Table 5: Overall SQIsign performance. KeyGen Cycle Speedup

Cycle

C ASM AVX-512 AVX-512 (Qlapoti)

81.54M 50.05M 46.21M

1.00× 1.63× 1.76×

185.29M 116.15M 108.05M

1.00× 1.60× 1.71×

12.54M 5.77M 3.94M

1.00× 2.17× 3.18×

28.11M

2.90×

68.81M

2.69×

3.94M

3.18×

C ASM AVX-512 AVX-512 (Qlapoti)

206.08M 152.80M 132.47M

1.00× 1.35× 1.56×

481.06M 349.64M 293.40M

1.00× 1.38× 1.64×

32.04M 20.39M 11.50M

1.00× 1.57× 2.79×

97.21M

2.12×

253.19M

1.90×

11.50M

2.79×

C ASM AVX-512 AVX-512 (Qlapoti)

355.81M 251.91M 197.17M

1.00× 1.41× 1.80×

820.28M 582.24M 466.21M

1.00× 1.41× 1.76×

64.19M 39.75M 21.63M

1.00× 1.61× 2.97×

151.48M

2.35×

370.87M

2.21×

21.63M

2.97×

Level

Impl.

I

III

V

Sign Speedup

Verify Cycle Speedup

Related documents

• Nora: Normalized Orthogonal Row Alignment for Scalable Matrix Optimizer#155302
• RODE: A Radial-Orthogonal Decoupled Engine for Optimization#662204
• RODE: A Radial-Orthogonal Decoupled Engine for Optimization#661154
• RODE: A Radial-Orthogonal Decoupled Engine for Optimization#480064
• Amortized Neural Optimization for Pre-Layout Signal Integrity Design Space Exploration using Differentiable Surrogates#266162