# Fast Implementation of Multiplication on Polynomial Rings

Authors: Boyu Wang; Haiying Gao; Fan Yang
Journal: Security and Communication Networks
 (2022)

---

## Abstract

Multiplication on polynomial rings has been widely used in public-key cryptographic schemes based on ideal lattices. It is an important module that significantly affects the efficiency of the schemes. Improved Preprocess-then-NTT (IPtNTT) is an algorithm which can fast realize multiplication on polynomial rings. Compared with the Number Theoretic Transform (NTT), the IPtNTT weakens the parameter restriction of lattice-based public-key cryptographic schemes. By optimizing the IPtNTT with the AVX2 instruction set, we reduce the clock cycles consumed by multiplication on polynomial rings to 15%–22%. According to the experimental results, we give specific suggestions on using AVX2 optimized IPtNTT to realize multiplication on polynomial rings with different parameters chosen in lattice-based public-key cryptosystems.


---

## Body

## 1. Introduction

Lattice-based public-key cryptography has emerged as a hot research direction in postquantum cryptography. Compared with the classic hard lattice problems such as the shortest vector problem and the closest vector problem, the Learning with Errors (LWE) problem [1] seems more suitable for constructing lattice-based public-key cryptographic schemes, especially the Ring-Learning with Errors (RLWE) problem [2] on which NEWHOPE is based, and the Module-Learning with Errors (MLWE) problem [3] on which KYBER is based. Multiplication on polynomial rings is a necessary and time-consuming operation in lattice cryptosystems. Through studying its fast implementation, we can significantly improve the multiplication efficiency.Various algorithms including Karatsuba algorithm [4, 5] and Fast Fourier Transform (FFT) algorithm [6] have been proposed to improve the implementation efficiency of polynomial multiplication. If we take the coefficients of the polynomial from the quotient ring ℤq=ℤ/qℤ, we call FFT the Number Theoretic Transform (NTT) [6]. NTT uses the root of unity in the residue class ring instead of that in the field of complexity number applied in FFT, which improves the calculation accuracy. In a polynomial ring Rq=ℤqx/xn+1, where q is a prime and n is some power of 2, the product of two polynomials corresponds to the cyclic convolution of the coefficient vectors. If q and n satisfy q≡1mod2n, which means there exist 2n-th roots of unity in ℤq, we can efficiently multiply polynomials through NTT.However, the condition thatq and n must satisfy to fasten the multiplication on polynomial rings through NTT severely limits the parameter selection in lattice-based public-key cryptosystem. In Inscrypt 2018, Zhou et al. [7] proposed a technique called Preprocess-then-NTT (PtNTT) which divides the n-th polynomials into several polynomials of lower degree n/2r-th, where r is a positive integer) and then performs the NTT for these n/2r-th polynomials, respectively. Since n and q are required to satisfy n/2r−1|q−1, the selection range of q in lattice-based public-key cryptosystem is enlarged. But the calculation efficiency of multiplication through PtNTT is lower than that through the original NTT due to some additional procedures.In order to solve the efficiency problem, Zhu et al. [8] introduced a technique called Improved Preprocess-then-NTT (IPtNTT). Similar to PtNTT, IPtNTT also divides the polynomials into several polynomials of lower degree and then applies NTT to their multiplications. But in the final procedure, they completed the multiplication using some ideas of Karatsuba algorithm, so as to optimize the implementation efficiency.To further improve the implementation efficiency, advanced vector extensions (AVX), a new set of SIMD instructions that are compatible with AMD’s SSE5 instruction sets, are often used. The AVX instructions deal with 128-bit wide data. Making better use of registers, the AVX2 instruction set widens the width to 256 bits to improve the computing capability of the processor. Seiler [9] presented a method to optimize NTT with AVX2 instructions to implement NTT more efficiently. However, how to use AVX2 instruction set to accelerate IPtNTT needs further research.Aiming at accelerating the implementation of IPtNTT algorithm with AVX2 instruction set, we first give the AVX2 implementation methods of modular operation and coefficients shuffling, respectively. We then demonstrate our AVX2 optimization toNTT and NTT−1 transformation on this basis. After deducing the specification of the reduction to the point-wise multiplications in [8], we finally present AVX2 optimized IPtNTT. Having made several experiments with some common parameters used in lattice-based public-key encryption schemes, we compare the implementation efficiency before and after optimization. The experimental results show that the clock cycles consumed by AVX2 optimized αIPtNTT are 15%–22% of standard C implementation. According to the experimental results, we give specific suggestions on using AVX2 optimized IPtNTT to realize multiplication on polynomial rings with certain q and n.


### 1.1. Roadmap

The remainder of the paper is organized as follows. In Section2, we provide preliminaries and notations. In Section 3, we present the AVX2 optimization to IPtNTT. In Section 4, we show the experimental results and the analysis. In Section 5, we give the conclusion.


## 1.1. Roadmap

The remainder of the paper is organized as follows. In Section2, we provide preliminaries and notations. In Section 3, we present the AVX2 optimization to IPtNTT. In Section 4, we show the experimental results and the analysis. In Section 5, we give the conclusion.


## 2. Preliminaries

Letℤ be the integer ring, q be a prime, and ℤq represent the residue class ring ℤ/qℤ, which means ℤq=0,…,q−1. We have Rq=ℤqx/xn+1 which denotes the polynomial ring consisting of polynomials modulo xn+1, whose coefficients are chosen from ℤq. For an arbitrary element f belonging to Rq, it can be written as f=∑i=0n−1fixi, where fi∈ℤq. We can also consider f as a vector containing coefficients of the polynomial, which means f↔fn−1,fn−2,…,f0. For f,g∈Rq, we have f∘g which denotes the point-wise multiplication of f and g, which means f∘g↔fn−1·gn−1modq,fn−2·gn−2modq,…,f0·g0modq. f·g represents the product or the cyclic convolution of f and g. If h=f·g∈Rq, then(1)hi=∑j=0ifjgi−j−∑j=i+1n−1fjgn+i−jmodq,and we can see its computational complexity is Tn=n2.


### 2.1. Multiplication on Polynomial Rings Based on NTT

We can fast realize the multiplication on polynomial rings through the forward transformationNTT, the inverse transformation NTT−1, and the point-wise multiplication on Rq [10]. We then present the corresponding basic knowledge.Now, we have an integern which satisfies n=2m, m∈ℤ, and a prime q which satisfies q≡1mod2n. Let ω be a n-th primitive root of unity in ℤq and γ be a 2n-th primitive root of unity, which means γ=ω. Then for ∀f∈Rq, we can describe its forward transformation NTT as f^i=∑j=0n−1γjfj·ωijmodq. Note that f^↔f^n−1,f^n−2,…,f^0, where f^=NTTf. Accordingly, for ∀f^∈Rq, the inverse transformation NTT−1 can be described as fi=n−1∑j=0n−1γjf^j·ω−ijmodq, which means f=NTT−1f^. So, we can compute f·g by(2)f·g=NTT−1NTTf°NTTg.If we precompute and store the powers of the root, the computational complexity ofNTT or NTT−1 transformation is Tn=n2, and that of equation (2) is Tn=3n2+n. It is an unacceptable complexity which is even larger than that of equation (1). However, Eleanor et al. [11] proposed a fast implementation to reduce the computational complexity of NTT and NTT−1 transformation to Tn=nlogn, thus reducing that of equation (2) to Tn=3nlogn+n.Here, we briefly introduce a method with divide-and-conquer strategy to implement NTT fast. Dividing the sequence of polynomial coefficients into twon/2 long subsequences according to the parity of their subscripts, we then have fi′=f2i and fi″=f2i+1, where i=0,1,…,n/2−1. Their NTT transformed results turn out to be f^i′=∑j=0n/2−1γ˜j·fj′·ω˜ijmodq and f^i″=∑j=0n/2−1γ˜j·fj″·ω˜ijmodq, respectively, where ω˜ is a n/2-th primitive root of unity in ℤq and γ˜ is a n-th root, which means γ˜=ω. After the transformation in this level, the computational complexity changes from one n-th NTT transformation to two n/2-th NTTs, which is Tn=n2/2. We repeat such transformation level by level until each subsequence contains only 1 coefficient and finally do the NTT to each pair of subsequences. In this way, the computational complexity is Tn=nlogn. Note that the subscripts of NTT transformed coefficients are bit-reversed. If k is an integer smaller than n which can be described as k=km−12m−1+km−22m−2+⋯+k12+k0 and brv k represents the bit-reversed k, we have brvk=k02m−1+k12m−2+⋯+km−22+km−1. In order to demonstrate the transformation more explicitly, we use fx to represent the x-th coefficient of f with x shown in binary while being transformed level by level. Let f_i be the transformation result in the i-th level and Abrvx be the brvx-th coefficient in the last level; then, the result of the m-th level is f^↔A0,A1,…,An−1↔f_mbrv0,f_mbrvn/2,…,f_mbrvn−1. Therefore, the transformation in the first level can be described as follows:(3)f_1i=f_0i+γbrv1f_0i+brv1,f_1i+brv1=f_0i−γbrv1f_0i+brv1,0≤i≤n2−1.Similarly, we can construct relationships of subsequent levels. Figure1 shows the structure of the fast NTT implementation. By inversely doing the transformation, we get the structure of NTT−1 which shares the same complexity with NTT.Figure 1
Structure of fastNTT implementation.


### 2.2. Improved Preprocess-Then-NTT

Algorithm 1: IPtNTT.
Input: coefficient vectors of f and gOutput: coefficient vector of h=f·gStep 1: (Separation)fx=∑i=02α−1xi·f¯ix2α, gx=∑j=02α−1xj·g¯jx2α, where f¯iz, g¯jz are all n/2α-th.Step 2: (Multiplication) Computeh¯iz∈ℤqz/zn/2α+1 as follows: h¯iz=∑l=0if¯lz·g¯i−lz+∑l=i+12α−1z·f¯lz·g¯2α+i−lz=NTT−1∑l=0iNTTf¯lz°NTTg¯i−lz+NTTz°∑l=i+12α−1NTTf¯lz°NTTg¯2α+i−lz, where 0≤i≤2α.Step 3: (Gatheration) Computeh in Rq as follows: hx=∑i=02α−1xi·h¯ix2α.In IPtNTT,f and g are divided into several shorter polynomials. After NTT, point-wise multiplication, and NTT−1, the coefficients of those shorter polynomials turn to be a series of coefficient vectors h¯i, which are finally combined into the coefficient vector corresponding to f·g. Algorithm 1 is the α-round IPtNTT (αIPtNTT) algorithm proposed by Zhu et al. [8], which requires q≡1modn/2α−1. In Algorithm 1, z represents x2α, where α≥1, and NTTz can be precomputed. In Step 2, when they compute each part, for i=j, they do it as usual, while for i≠j, they compute two corresponding parts together as follows in practice:(4)NTTf¯i°NTTg¯j+NTTf¯j°NTTg¯i=NTTf¯i+NTTf¯j°NTTg¯i+NTTg¯j−NTTf¯i°NTTg¯i−NTTf¯j°NTTg¯j.In this way, they reduce the point-wise multiplications wheni≠j. Throughout this paper, we use ∗^ to represent NTT∗.Their result shows thatαIPtNTT is composed of n/2α-th NTTs, n/2α-th NTT−1s, and n/2α-th point-wise multiplications, and the number of them is 2α+1, 2α and 3·22α−2+2α−1, respectively. So, the total computational complexity is as follows:(5)Tn=3nlogn+3·2α−2−3α+12n.


## 2.1. Multiplication on Polynomial Rings Based on NTT

We can fast realize the multiplication on polynomial rings through the forward transformationNTT, the inverse transformation NTT−1, and the point-wise multiplication on Rq [10]. We then present the corresponding basic knowledge.Now, we have an integern which satisfies n=2m, m∈ℤ, and a prime q which satisfies q≡1mod2n. Let ω be a n-th primitive root of unity in ℤq and γ be a 2n-th primitive root of unity, which means γ=ω. Then for ∀f∈Rq, we can describe its forward transformation NTT as f^i=∑j=0n−1γjfj·ωijmodq. Note that f^↔f^n−1,f^n−2,…,f^0, where f^=NTTf. Accordingly, for ∀f^∈Rq, the inverse transformation NTT−1 can be described as fi=n−1∑j=0n−1γjf^j·ω−ijmodq, which means f=NTT−1f^. So, we can compute f·g by(2)f·g=NTT−1NTTf°NTTg.If we precompute and store the powers of the root, the computational complexity ofNTT or NTT−1 transformation is Tn=n2, and that of equation (2) is Tn=3n2+n. It is an unacceptable complexity which is even larger than that of equation (1). However, Eleanor et al. [11] proposed a fast implementation to reduce the computational complexity of NTT and NTT−1 transformation to Tn=nlogn, thus reducing that of equation (2) to Tn=3nlogn+n.Here, we briefly introduce a method with divide-and-conquer strategy to implement NTT fast. Dividing the sequence of polynomial coefficients into twon/2 long subsequences according to the parity of their subscripts, we then have fi′=f2i and fi″=f2i+1, where i=0,1,…,n/2−1. Their NTT transformed results turn out to be f^i′=∑j=0n/2−1γ˜j·fj′·ω˜ijmodq and f^i″=∑j=0n/2−1γ˜j·fj″·ω˜ijmodq, respectively, where ω˜ is a n/2-th primitive root of unity in ℤq and γ˜ is a n-th root, which means γ˜=ω. After the transformation in this level, the computational complexity changes from one n-th NTT transformation to two n/2-th NTTs, which is Tn=n2/2. We repeat such transformation level by level until each subsequence contains only 1 coefficient and finally do the NTT to each pair of subsequences. In this way, the computational complexity is Tn=nlogn. Note that the subscripts of NTT transformed coefficients are bit-reversed. If k is an integer smaller than n which can be described as k=km−12m−1+km−22m−2+⋯+k12+k0 and brv k represents the bit-reversed k, we have brvk=k02m−1+k12m−2+⋯+km−22+km−1. In order to demonstrate the transformation more explicitly, we use fx to represent the x-th coefficient of f with x shown in binary while being transformed level by level. Let f_i be the transformation result in the i-th level and Abrvx be the brvx-th coefficient in the last level; then, the result of the m-th level is f^↔A0,A1,…,An−1↔f_mbrv0,f_mbrvn/2,…,f_mbrvn−1. Therefore, the transformation in the first level can be described as follows:(3)f_1i=f_0i+γbrv1f_0i+brv1,f_1i+brv1=f_0i−γbrv1f_0i+brv1,0≤i≤n2−1.Similarly, we can construct relationships of subsequent levels. Figure1 shows the structure of the fast NTT implementation. By inversely doing the transformation, we get the structure of NTT−1 which shares the same complexity with NTT.Figure 1
Structure of fastNTT implementation.


## 2.2. Improved Preprocess-Then-NTT

Algorithm 1: IPtNTT.
Input: coefficient vectors of f and gOutput: coefficient vector of h=f·gStep 1: (Separation)fx=∑i=02α−1xi·f¯ix2α, gx=∑j=02α−1xj·g¯jx2α, where f¯iz, g¯jz are all n/2α-th.Step 2: (Multiplication) Computeh¯iz∈ℤqz/zn/2α+1 as follows: h¯iz=∑l=0if¯lz·g¯i−lz+∑l=i+12α−1z·f¯lz·g¯2α+i−lz=NTT−1∑l=0iNTTf¯lz°NTTg¯i−lz+NTTz°∑l=i+12α−1NTTf¯lz°NTTg¯2α+i−lz, where 0≤i≤2α.Step 3: (Gatheration) Computeh in Rq as follows: hx=∑i=02α−1xi·h¯ix2α.In IPtNTT,f and g are divided into several shorter polynomials. After NTT, point-wise multiplication, and NTT−1, the coefficients of those shorter polynomials turn to be a series of coefficient vectors h¯i, which are finally combined into the coefficient vector corresponding to f·g. Algorithm 1 is the α-round IPtNTT (αIPtNTT) algorithm proposed by Zhu et al. [8], which requires q≡1modn/2α−1. In Algorithm 1, z represents x2α, where α≥1, and NTTz can be precomputed. In Step 2, when they compute each part, for i=j, they do it as usual, while for i≠j, they compute two corresponding parts together as follows in practice:(4)NTTf¯i°NTTg¯j+NTTf¯j°NTTg¯i=NTTf¯i+NTTf¯j°NTTg¯i+NTTg¯j−NTTf¯i°NTTg¯i−NTTf¯j°NTTg¯j.In this way, they reduce the point-wise multiplications wheni≠j. Throughout this paper, we use ∗^ to represent NTT∗.Their result shows thatαIPtNTT is composed of n/2α-th NTTs, n/2α-th NTT−1s, and n/2α-th point-wise multiplications, and the number of them is 2α+1, 2α and 3·22α−2+2α−1, respectively. So, the total computational complexity is as follows:(5)Tn=3nlogn+3·2α−2−3α+12n.


## 3. AVX2 Optimization to IPtNTT

Taking advantage of its parallel computation, we use AVX2 instructions to optimize IPtNTT. As the parameterqs commonly chosen by LWE cryptosystems like 3329, 7681, 12 289 or the even smaller one 1409 are all within 16 bits, we determine to store 16 numbers in 1 register. So, each AVX2 instruction deals with vectors containing 16 numbers. According to Algorithm 1, we find that the key events to optimize αIPtNTT include AVX2 optimized modular reduction, coefficients shuffling, NTT, NTT−1, and point-wise multiplication. The logical relationship is shown in Figure 2. We introduce our implementation, respectively, in the following.Figure 2
Modules of AVX2 optimized IPtNTT.


### 3.1. AVX2 Implemented Modular Reduction

As there exists no instruction to acquire remainders, we implement the modified Montgomery transformation and the modified Barrett transformation proposed by Seiler [9] which are called Montgomery transformation and Barrett transformation, respectively, throughout this paper. Specifically, Montgomery transformation deals with numbers within 32 bits, while Barrett transformation deals with those within 16 bits. This subsection presents how we implement such transformations with AVX2 instructions.As for the specific realization of the algorithm, we move the product of two 16-bit numbers to the right for 16 bits by using the instruction which gets the high 16 bits of the product. The instruction to acquire low 16 bits can efficiently get the remainder of the product modulo216. In Barrett transformation, a number that can be precomputed needs to be rounded, which can be achieved by adding 0.5 to the original number and then rounding it down. The specific implementation procedure is shown in Algorithms 2 and 3.Algorithm 2:Montgomery transformation.
Input: An integer a satisfies −215q<a=216a1+a0<215q, where 0≤a0<216, q<215Output: r′=2−16amodq, where −q<r′<qStep 1: Using _mm256_mullo_epi16 (∗,∗∗), calculate and store the lower 16 bits of the product of two numbers;Step 2: Using _mm256_mulhi_epi16 (∗,∗∗), calculate and store the higher 16 bits of the signed product of two numbers;Step 3: Using _mm256_mullo_epi16 (∗, q−1), compute and store the lower 16 bits of the product between the result of Step 1 and q−1, which has been precomputed;Step 4: Using _mm256_mulhi_epi16 (∗, q), calculate and store the higher 16 bits of the product of the result of Step 3 and q;Step 5: Using _mm256_sub_epi16 (∗,∗∗), calculate and store the difference between the result of Step 2 and the result of Step 4.Algorithm 3:Barrett transformation.
Input: An integer a satisfies −215<a<215, where 0≤q<215Output:r=amodq, where 0≤r≤qStep 1: Using _mm256_mulhi_epi16 (a, v), calculate and store the high 16 bits of the signed product of a and v, where v, which is the integer closest to 2⌊logq⌋+15/q, has been precomputed;Step 2: Using _mm256_srai_epi16 (∗, ⌊logq⌋−1), shift the result of Step 1 to the right for ⌊logq⌋−1 bits;Step 3: Using _mm256_mullo_epi16 (∗, q), calculate and store the lower 16 bits of the product of the result of Step 2 and q;Step 4: Using _mm256_sub_epi16 (a,∗), calculate and store the difference between a and the result of Step 3.


### 3.2. AVX2 Implemented Shuffling

According to the description in 2.1, the fast implementation of NTT is hierarchical. In the forward transformation NTT, for each level before the fourth level from the bottom, all sixteen coefficients in one vector use the same power of the chosen root. So, we can use _mm256_loadu_si256 to take out the two coefficient vectors and directly do the computation according to the relationship of each two levels similar to equation (3). In the last 4 levels, after taken out by _mm256_loadu_si256, coefficients in the two vectors need different powers of the root to finish the transformation, but the two coefficients in one relation stored in the same vector. Therefore, it is necessary to shuffle the elements of the coefficient vectors.In the fourth level from the bottom in NTT transformation, eight coefficients in one vector are performed with the other eight using the same root. We need to switch the lower eight in the vector containing higher subscript coefficients with the higher eight in the lower vector. We first use _mm256_permute4x64_epi64 (∗,0xD8) to both vectors and then get the coefficients with higher subscripts in those relations with _mm256_unpackhi_epi64 (∗,∗∗) and the lowers with _mm256_unpacklo_epi64 (∗,∗∗). After we use _mm256_permute4x64_epi64 (∗,0xD8) to both vectors again, the two corresponding coefficients in relations similar to equation (3) are at the corresponding positions of the two vectors. We present the process of coefficients shuffling in the fourth level from the bottom of the forward transformation in Figure 3. Notice that there is no need to restore the coefficients to the original order for the next level after NTT. We can do above shuffles the other way around in the first four levels in NTT−1. Shuffles in other levels can be performed similarly.Figure 3
Shuffle in the fourth level from the bottom in NTT, where scr1 and scr2 are two vectors.Algorithm 4:Implementation of NTT.
Input: Coefficient vector of f, which is f0,f1,…,fn−1Output: Bit-reversed NTT coefficient vector of f, which is f^brv0,f^brv1,…,f^brvn−1Step 1: Using _mm256_set1_epi16 (∗), or _mm256_set_epi16 for the last 4 levels, store the value from the powers table of the root used in this round;Step 2: Using _mm256_load_si256 (∗), store the coefficient vectors;Step 3: Shuffle if it is in the last 4 levels, skip this otherwise;Step 4: Montgomery transform the products of higher degree coefficients and the result of Step 1;Step 5: Using _mm256_sub_epi16 (∗,∗∗), compute the difference between the lower degree coefficients and the result of Step 4 to the higher degree coefficient vector;Step 6: Using _mm256_add_epi16 (∗,∗∗), compute the sum of the lower degree coefficients and the result of Step 4 to the lower degree coefficient vector;Step 7: Barrett transform the results of Step 5 and Step 6, respectively;Step 8: Restore the coefficients if it is in the last 4 levels, skip this otherwise;Step 9: Using _mm256_store_si256 (∗), store the result of the last step;Step 10: Circulate from Step 2, until the value stored in Step 1 needs to be updated;Step 11: Circulate from Step 1, until all coefficients have been transformed;Step 12: Enter the next level, repeat steps above till the last level.


### 3.3. AVX2 Implementation of NTT andNTT−1

As mentioned above, the fast implementation of NTT relies on hierarchical operations, each of which involves coefficients’ modulo multiplication with the corresponding powers of the root of unity. Since we use Montgomery transformation in our AVX2 implementation, what we actually get is2−16amodq, where a is our target. Therefore, when we build the table of the powers of the root, we multiply the brvi-th power of the root with 216 and store the remainder of the product modulo q. Then, the Montgomery transformed product of any integer and the value in the table is within (−q, q). The powers table of the inverse of the root used for the inverse transformation can also be constructed in this way. To ensure the correctness of the results, we additionally perform Barrett transformation on the results of each level. The procedures of the forward transformation NTT and the inverse transformation NTT−1 are presented in Algorithms 4 and 5, respectively.When doing the point-wise multiplication, we only need to use AVX2 instructions to do the parallel computation. It should be noted that those multiplications which cannot be precomputed can be implemented through two Montgomery transformations. To achieve this, first Montgomery transform the product of a multiplier and232modq, and then Montgomery transform the product of the result and the other multiplier.Algorithm 5:Implementation of NTT−1.
Input: Bit-reversed NTT coefficient vector of f, which is f^brv0,f^brv1,…,f^brvn−1Output: Coefficient vector of f, which is f0,f1,…,fn−1Step 1: Using _mm256_set1_epi16 (∗), or _mm256_set_epi16 for the first 4 levels, store the value from the powers table of the root used in this round;Step 2: Using _mm256_load_si256 (∗), store the coefficient vectors;Step 3: Shuffle if it is in the last 4 levels, skip this otherwise;Step 4: Using _mm256_add_epi16 (∗,∗∗), compute the sum of the two coefficient vectors from Step 3 or Step 2 to the lower degree coefficient vector;Step 5: Using _mm256_sub_epi16 (∗,∗∗), compute the difference of the lower degree coefficient vector and the higher one from Step 3 or Step 2;Step 6: Montgomery transform the products of the higher degree coefficients and the result of Step 1 and assign them to the higher degree coefficient vector;Step 7: Barrett transform the result of Step 4 and Step 6, respectively;Step 8: Restore the coefficients if it is in the first 4 levels, skip this otherwise;Step 9: Using _mm256_store_si256 (∗), store the result of the last step;Step 10: Circulate from Step 2, until the value stored in Step 1 needs to be updated;Step 11: Circulate from Step 1, until all coefficients have been transformed;Step 12: Enter the next level, repeat steps above till the last level.Step 13: Montgomery transform the products of the result of Step 12 and216n−1modq, then do the Barrett transformation.


### 3.4. AVX2 Implementation of IPtNTT

In order to improve the computation efficiency, we take advantages of some techniques in point-wise multiplication. The technical details of IPtNTT are given below.From 2.2, we can see the time complexity of IPtNTT mainly includes two parts. One is the complexity of NTT andNTT−1 to the shorter polynomials, i.e., f¯iz and g¯iz, the AVX2 implementation of which has been presented in 3.3. The other is the complexity of the point-wise multiplications mentioned in Step 2 of Algorithm 1. We deduce an applicable expression for arbitrary rounds of IPtNTT based on [8]. In the situation that we have acquired f¯^i and g¯^i, i.e., NTT transformed coefficients of f¯iz and g¯iz, where 0≤i<2α, we may acquire h¯^i, the NTT coefficients of h¯i as follows:(6)i≡0mod2:h¯^i=f¯^i2°g¯^i2+∑l=0i2−1f¯^l+f¯^i−l°g¯^l+g¯^i−l−f¯^l°g¯^l−f¯^i−l°g¯^i−l+z^°f¯^2α2+i2°g¯^2α2+i2+∑l=i+12α+i2−1f¯^l+f¯^2α+i−l°g¯^l+g¯^2α+i−l−f¯^l°g¯^l−f¯^2α+i−l°g¯^2α+i−li≡1mod2:h¯^i=∑l=0i−12f¯^l+f¯^i−l°g¯^l+g¯^i−l−f¯^l°g¯^l−f¯^i−l°g¯^i−l+z^°∑l=i+12α+i−12f¯^l+f¯^2α+i−l°g¯^l+g¯^2α+i−l−f¯^l°g¯^l−f¯^2α+i−l°g¯^2α+i−l,It is clear thatf¯^i°g¯^i and f¯^i+f¯^j°g¯^i+g¯^j, where 0≤i,j<2α, can be used over and over again once computed to acquire h¯^i, so as to reduce the point-wise multiplications. Since the whole process includes 2α multiplications like f¯^i°g¯^i, C2α2 multiplications like f¯^i+f¯^j°g¯^i+g¯^j, and 2α−1 multiplications like z^°∗, αIPtNTT includes 22α−1+3·2α−1−1 point-wise multiplications in total. Including 2α+1NTTs and 2αNTT−1s, the computational complexity of αIPtNTT is as follows:(7)Tn=22α−1+3·2α−1−1·n2α+2α+1·n2αlogn2α+2α·n2αlogn2α=2α−1n+32n−n2α−3αn+3nlogn.Letα traverse ℤ10, and we compute Tn when n=256,512,1024, and present the result in Table 1. Note that z^ can also be precomputed as z^=NTT0,1,0,…,0, and we store the remainder of the product of z^ and 216 modulo q for the convenience of Montgomery transformation. The AVX2 implementation equation (6) can be described as Algorithm 6.Table 1
Tn of IPtNTT with different n and q.
αn = 256n = 512n = 10240640014 33631 7441588813 31229 6962544012 41627 9043521611 96827 0084548812 51228 0965677615 08833 248610 10821 75246 576717 53436 60476 280833 15167 838138 748965 151.5131 839266 750Algorithm 6:Implementation of point-wise multiplication in IPtNTT.
Input: f¯^i and g¯^i, NTT coefficient vectors of f¯iz and g¯iz, where 0≤i<2αOutput: h¯^i, NTT coefficient vector of h¯iStep 1: Montgomery transform the product off¯^i and 232modq;Step 2: Montgomery transform the result off¯^i°g¯^i and f¯^i+f¯^j°g¯^i+g¯^j, where f¯^i is replaced by the result of Step 1;Step 3: Compute and Barrett transform the part in parentheses which will perform point-wise multiplication withz^ in equation (6);Step 4: Montgomery transformz^°∗ in equation (6);Step 5: Compute and Barrett transform the remaining additions and subtractions in equation (6).Considering the limitation to the input, we need to do different times of Barrett transformations as the chosenq differs. For example, if q=12289, we need to Barrett transform after every addition or subtraction; if q=7681, we need to Barrett transform after every three additions or subtractions; if q=3329 or even less, we do not need to do extra Barrett transformations when α≤3.


## 3.1. AVX2 Implemented Modular Reduction

As there exists no instruction to acquire remainders, we implement the modified Montgomery transformation and the modified Barrett transformation proposed by Seiler [9] which are called Montgomery transformation and Barrett transformation, respectively, throughout this paper. Specifically, Montgomery transformation deals with numbers within 32 bits, while Barrett transformation deals with those within 16 bits. This subsection presents how we implement such transformations with AVX2 instructions.As for the specific realization of the algorithm, we move the product of two 16-bit numbers to the right for 16 bits by using the instruction which gets the high 16 bits of the product. The instruction to acquire low 16 bits can efficiently get the remainder of the product modulo216. In Barrett transformation, a number that can be precomputed needs to be rounded, which can be achieved by adding 0.5 to the original number and then rounding it down. The specific implementation procedure is shown in Algorithms 2 and 3.Algorithm 2:Montgomery transformation.
Input: An integer a satisfies −215q<a=216a1+a0<215q, where 0≤a0<216, q<215Output: r′=2−16amodq, where −q<r′<qStep 1: Using _mm256_mullo_epi16 (∗,∗∗), calculate and store the lower 16 bits of the product of two numbers;Step 2: Using _mm256_mulhi_epi16 (∗,∗∗), calculate and store the higher 16 bits of the signed product of two numbers;Step 3: Using _mm256_mullo_epi16 (∗, q−1), compute and store the lower 16 bits of the product between the result of Step 1 and q−1, which has been precomputed;Step 4: Using _mm256_mulhi_epi16 (∗, q), calculate and store the higher 16 bits of the product of the result of Step 3 and q;Step 5: Using _mm256_sub_epi16 (∗,∗∗), calculate and store the difference between the result of Step 2 and the result of Step 4.Algorithm 3:Barrett transformation.
Input: An integer a satisfies −215<a<215, where 0≤q<215Output:r=amodq, where 0≤r≤qStep 1: Using _mm256_mulhi_epi16 (a, v), calculate and store the high 16 bits of the signed product of a and v, where v, which is the integer closest to 2⌊logq⌋+15/q, has been precomputed;Step 2: Using _mm256_srai_epi16 (∗, ⌊logq⌋−1), shift the result of Step 1 to the right for ⌊logq⌋−1 bits;Step 3: Using _mm256_mullo_epi16 (∗, q), calculate and store the lower 16 bits of the product of the result of Step 2 and q;Step 4: Using _mm256_sub_epi16 (a,∗), calculate and store the difference between a and the result of Step 3.


## 3.2. AVX2 Implemented Shuffling

According to the description in 2.1, the fast implementation of NTT is hierarchical. In the forward transformation NTT, for each level before the fourth level from the bottom, all sixteen coefficients in one vector use the same power of the chosen root. So, we can use _mm256_loadu_si256 to take out the two coefficient vectors and directly do the computation according to the relationship of each two levels similar to equation (3). In the last 4 levels, after taken out by _mm256_loadu_si256, coefficients in the two vectors need different powers of the root to finish the transformation, but the two coefficients in one relation stored in the same vector. Therefore, it is necessary to shuffle the elements of the coefficient vectors.In the fourth level from the bottom in NTT transformation, eight coefficients in one vector are performed with the other eight using the same root. We need to switch the lower eight in the vector containing higher subscript coefficients with the higher eight in the lower vector. We first use _mm256_permute4x64_epi64 (∗,0xD8) to both vectors and then get the coefficients with higher subscripts in those relations with _mm256_unpackhi_epi64 (∗,∗∗) and the lowers with _mm256_unpacklo_epi64 (∗,∗∗). After we use _mm256_permute4x64_epi64 (∗,0xD8) to both vectors again, the two corresponding coefficients in relations similar to equation (3) are at the corresponding positions of the two vectors. We present the process of coefficients shuffling in the fourth level from the bottom of the forward transformation in Figure 3. Notice that there is no need to restore the coefficients to the original order for the next level after NTT. We can do above shuffles the other way around in the first four levels in NTT−1. Shuffles in other levels can be performed similarly.Figure 3
Shuffle in the fourth level from the bottom in NTT, where scr1 and scr2 are two vectors.Algorithm 4:Implementation of NTT.
Input: Coefficient vector of f, which is f0,f1,…,fn−1Output: Bit-reversed NTT coefficient vector of f, which is f^brv0,f^brv1,…,f^brvn−1Step 1: Using _mm256_set1_epi16 (∗), or _mm256_set_epi16 for the last 4 levels, store the value from the powers table of the root used in this round;Step 2: Using _mm256_load_si256 (∗), store the coefficient vectors;Step 3: Shuffle if it is in the last 4 levels, skip this otherwise;Step 4: Montgomery transform the products of higher degree coefficients and the result of Step 1;Step 5: Using _mm256_sub_epi16 (∗,∗∗), compute the difference between the lower degree coefficients and the result of Step 4 to the higher degree coefficient vector;Step 6: Using _mm256_add_epi16 (∗,∗∗), compute the sum of the lower degree coefficients and the result of Step 4 to the lower degree coefficient vector;Step 7: Barrett transform the results of Step 5 and Step 6, respectively;Step 8: Restore the coefficients if it is in the last 4 levels, skip this otherwise;Step 9: Using _mm256_store_si256 (∗), store the result of the last step;Step 10: Circulate from Step 2, until the value stored in Step 1 needs to be updated;Step 11: Circulate from Step 1, until all coefficients have been transformed;Step 12: Enter the next level, repeat steps above till the last level.


## 3.3. AVX2 Implementation of NTT andNTT−1

As mentioned above, the fast implementation of NTT relies on hierarchical operations, each of which involves coefficients’ modulo multiplication with the corresponding powers of the root of unity. Since we use Montgomery transformation in our AVX2 implementation, what we actually get is2−16amodq, where a is our target. Therefore, when we build the table of the powers of the root, we multiply the brvi-th power of the root with 216 and store the remainder of the product modulo q. Then, the Montgomery transformed product of any integer and the value in the table is within (−q, q). The powers table of the inverse of the root used for the inverse transformation can also be constructed in this way. To ensure the correctness of the results, we additionally perform Barrett transformation on the results of each level. The procedures of the forward transformation NTT and the inverse transformation NTT−1 are presented in Algorithms 4 and 5, respectively.When doing the point-wise multiplication, we only need to use AVX2 instructions to do the parallel computation. It should be noted that those multiplications which cannot be precomputed can be implemented through two Montgomery transformations. To achieve this, first Montgomery transform the product of a multiplier and232modq, and then Montgomery transform the product of the result and the other multiplier.Algorithm 5:Implementation of NTT−1.
Input: Bit-reversed NTT coefficient vector of f, which is f^brv0,f^brv1,…,f^brvn−1Output: Coefficient vector of f, which is f0,f1,…,fn−1Step 1: Using _mm256_set1_epi16 (∗), or _mm256_set_epi16 for the first 4 levels, store the value from the powers table of the root used in this round;Step 2: Using _mm256_load_si256 (∗), store the coefficient vectors;Step 3: Shuffle if it is in the last 4 levels, skip this otherwise;Step 4: Using _mm256_add_epi16 (∗,∗∗), compute the sum of the two coefficient vectors from Step 3 or Step 2 to the lower degree coefficient vector;Step 5: Using _mm256_sub_epi16 (∗,∗∗), compute the difference of the lower degree coefficient vector and the higher one from Step 3 or Step 2;Step 6: Montgomery transform the products of the higher degree coefficients and the result of Step 1 and assign them to the higher degree coefficient vector;Step 7: Barrett transform the result of Step 4 and Step 6, respectively;Step 8: Restore the coefficients if it is in the first 4 levels, skip this otherwise;Step 9: Using _mm256_store_si256 (∗), store the result of the last step;Step 10: Circulate from Step 2, until the value stored in Step 1 needs to be updated;Step 11: Circulate from Step 1, until all coefficients have been transformed;Step 12: Enter the next level, repeat steps above till the last level.Step 13: Montgomery transform the products of the result of Step 12 and216n−1modq, then do the Barrett transformation.


## 3.4. AVX2 Implementation of IPtNTT

In order to improve the computation efficiency, we take advantages of some techniques in point-wise multiplication. The technical details of IPtNTT are given below.From 2.2, we can see the time complexity of IPtNTT mainly includes two parts. One is the complexity of NTT andNTT−1 to the shorter polynomials, i.e., f¯iz and g¯iz, the AVX2 implementation of which has been presented in 3.3. The other is the complexity of the point-wise multiplications mentioned in Step 2 of Algorithm 1. We deduce an applicable expression for arbitrary rounds of IPtNTT based on [8]. In the situation that we have acquired f¯^i and g¯^i, i.e., NTT transformed coefficients of f¯iz and g¯iz, where 0≤i<2α, we may acquire h¯^i, the NTT coefficients of h¯i as follows:(6)i≡0mod2:h¯^i=f¯^i2°g¯^i2+∑l=0i2−1f¯^l+f¯^i−l°g¯^l+g¯^i−l−f¯^l°g¯^l−f¯^i−l°g¯^i−l+z^°f¯^2α2+i2°g¯^2α2+i2+∑l=i+12α+i2−1f¯^l+f¯^2α+i−l°g¯^l+g¯^2α+i−l−f¯^l°g¯^l−f¯^2α+i−l°g¯^2α+i−li≡1mod2:h¯^i=∑l=0i−12f¯^l+f¯^i−l°g¯^l+g¯^i−l−f¯^l°g¯^l−f¯^i−l°g¯^i−l+z^°∑l=i+12α+i−12f¯^l+f¯^2α+i−l°g¯^l+g¯^2α+i−l−f¯^l°g¯^l−f¯^2α+i−l°g¯^2α+i−l,It is clear thatf¯^i°g¯^i and f¯^i+f¯^j°g¯^i+g¯^j, where 0≤i,j<2α, can be used over and over again once computed to acquire h¯^i, so as to reduce the point-wise multiplications. Since the whole process includes 2α multiplications like f¯^i°g¯^i, C2α2 multiplications like f¯^i+f¯^j°g¯^i+g¯^j, and 2α−1 multiplications like z^°∗, αIPtNTT includes 22α−1+3·2α−1−1 point-wise multiplications in total. Including 2α+1NTTs and 2αNTT−1s, the computational complexity of αIPtNTT is as follows:(7)Tn=22α−1+3·2α−1−1·n2α+2α+1·n2αlogn2α+2α·n2αlogn2α=2α−1n+32n−n2α−3αn+3nlogn.Letα traverse ℤ10, and we compute Tn when n=256,512,1024, and present the result in Table 1. Note that z^ can also be precomputed as z^=NTT0,1,0,…,0, and we store the remainder of the product of z^ and 216 modulo q for the convenience of Montgomery transformation. The AVX2 implementation equation (6) can be described as Algorithm 6.Table 1
Tn of IPtNTT with different n and q.
αn = 256n = 512n = 10240640014 33631 7441588813 31229 6962544012 41627 9043521611 96827 0084548812 51228 0965677615 08833 248610 10821 75246 576717 53436 60476 280833 15167 838138 748965 151.5131 839266 750Algorithm 6:Implementation of point-wise multiplication in IPtNTT.
Input: f¯^i and g¯^i, NTT coefficient vectors of f¯iz and g¯iz, where 0≤i<2αOutput: h¯^i, NTT coefficient vector of h¯iStep 1: Montgomery transform the product off¯^i and 232modq;Step 2: Montgomery transform the result off¯^i°g¯^i and f¯^i+f¯^j°g¯^i+g¯^j, where f¯^i is replaced by the result of Step 1;Step 3: Compute and Barrett transform the part in parentheses which will perform point-wise multiplication withz^ in equation (6);Step 4: Montgomery transformz^°∗ in equation (6);Step 5: Compute and Barrett transform the remaining additions and subtractions in equation (6).Considering the limitation to the input, we need to do different times of Barrett transformations as the chosenq differs. For example, if q=12289, we need to Barrett transform after every addition or subtraction; if q=7681, we need to Barrett transform after every three additions or subtractions; if q=3329 or even less, we do not need to do extra Barrett transformations when α≤3.


## 4. Experimental Results and Analysis

We obtained the results of computational cost experiments of C implemented multiplication on polynomial rings by using NTT, 1IPtNTT, 2IPtNTT, and 3IPtNTT, which are marked as NTT_mul andαIPtNTT_mul, and their AVX2 optimization, which are marked as AVX2_NTT and AVX2_αIPtNTT, by running the implementation on a 2.5 GHZ Intel Core i5 processor. We compiled the C implementation and AVX2 optimization with gcc 7.5.0 on Ubuntu 18.04 with the core Linux 4.15.0. We chose several pairs of (n, q) in our experiments. For each pair, we performed 10 rounds of experiments with 216 times of tests in each round and compute the average cost of 1 test. We chose 1409, 3329, 7681, and 12 289 for q and 256, 512, and 1024 for n, respectively. We present our experimental results in Table 2.Table 2
Clock cycle cost of NTT, 1PtNTT, 2PtNTT, and 3PtNTT, where NS indicates that the parameters pair does not satisfy the condition to perform the algorithm.
q14093329768112289140933297681122893329768112289n256256256256512512512512102410241024NTT_mulNSNS23031.406605423016.4419296NSNSNS46541.6980948NSNS95212.301562AVX_NTTNSNS3469.33292483476.2602384NSNSNS7149.018152NSNS14670.0246522AVX_NTT/NTT_mulNSNS15.063%15.103%NSNSNS15.360%NSNS15.408%1IPtNTTNS21286.93662120217.803420620163.5064042NSNS42091.40263441854.439682NSNS84812.5423236AVX_1IPtNTTNS3538.69379663558.5259023507.3187116NSNS7148.15605247214.228739NSNS14728.800009AVX_1IPt/1PtNTTNS16.624%17.601%17.394%NSNS16.982%17.236%NSNS17.366%2IPtNTT20162.846478421023.092182819855.843655819957.357791NS43756.738741440770.929254240781.1600114NS85377.21374784782.4298712AVX_2IPtNTT3687.39071923670.39667623708.23750823842.9119612NS7183.07490087166.1647797353.4273076NS14754.159469814901.6090182AVX_2IPt/2IPtNTT18.288%17.459%18.676%19.256%NS16.416%17.577%18.031%NS17.281%17.576%3IPtNTT20804.729220421383.156091420668.62810720753.167272241250.49313743162.534684840712.61546440770.658417689402.829789884598.859875483725.4423572AVX_3IPtNTT3313.61332363415.15797563464.3434283677.48801767925.8995887941.37160568054.65378188680.374063615533.764530215840.095983216816.2292356AVX_3IPt/3IPtNTT15.927%15.971%16.761%17.720%19.214%18.399%19.784%21.291%17.375%18.724%20.085%Analyzing Table2, we may confirm that the computational cost of NTT, 1IPtNTT, 2IPtNTT, and 3IPtNTT after AVX2 optimization is reduced to 15%–22%. Thus, for certain q and n, when n = 512, 1024, we suggest implementing AVX_αIPtNTT, where α=minr:q≡1modn/2r−1,r≥0. For example, if q = 3329, n = 512, we suggest using AVX_2IPtNTT; if q = 7681, n = 512, we suggest using AVX_1IPtNTT; if q  = 12 289, n = 512, 1024, we suggest using AVX_NTT. When n = 256, AVX_3IPtNTT is the best choice if q is within 13 bits; if q = 12 289, AVX_NTT is the most efficient.


## 5. Conclusion

IPtNTT is a method to fast realize multiplication on polynomial rings. Due to its significant advantage in expanding the parameter selection range of lattice public-key cryptographic schemes, it can be widely used in the programming of lattice public-key cryptographic schemes. In this paper, we present the AVX2 optimization implementation method of IPtNTT and show the experimental results with different parameters, based on the AVX2 implementation of technique modules such as modular reduction, coefficients shuffling, and point-wise multiplication of coefficient vectors. According to the experimental results, the implementation suggestions on using IPtNTT to complete the multiplication on polynomial rings are given. The research results in this paper can be applied to the fast implementation of lattice-based public-key cryptography schemes.



---

*Source: 4649158-2022-02-14.xml*