Jean-Marc Robert
email: jean-marc.robert@univ-tln.fr

Pascal Véron
email: veron@univ-tln.fr

J.-M Robert

Faster Multiplication over F 2 [X] using AVX512 instruction set and VPCLMULQDQ instruction

Keywords: Finite field multiplication, Karatsuba, Toom-Cook, post-quantum cryptography, code based cryptography, AVX2, AVX512, VPCLMULQDQ

. In this work, we revisit the different existing constant-time algorithms for arbitrary polynomial multiplication. We explore the different Karatsuba and Toom-Cook constructions in order to determine the best combinations for each polynomial degree range, in the context of AVX2 and AVX512 instruction sets. This leads to different kernels and constructions in each case. In particular, in the context of AVX512, we use the VPCLMULQDQ instruction, which is a vectorized binary polynomial multiplication instruction. This instruction deals with up to four polynomial (of degree up to 63) multiplications, that is four operand pairs of 64-bit words with 128-bit word storing each results, the four results being stored in one single 512-bit word. This allows to divide by roughly 3 the retired instruction number of the operation in comparison with the AVX2 instruction set implementations, while the speedup is

Introduction

In 2017, the NIST launched a consultation dealing with the so-called "Post-Quantum Cryptography" (PQC) [START_REF]Post-Quantum Cryptography[END_REF], leading to think that a practical quantum computer might appear in the next two or three decades. Among the candidates known to resist against quantum computers, several submissions on the NIST platform are code based protocols. The public key cryptosystem of Mc Eliece marked the beginning of code based cryptography [START_REF] Mceliece | A Public-Key Cryptosystem Based On Algebraic Coding Theory[END_REF]. The security of most of the code based protocols relies on a decision problem which can be stated without using the terminology of coding theory: the SD (Syndrome Decoding) problem.

Input

: H a (k, n) matrix over F 2 , s ∈ F k 2 a column vector, p an integer. Question : Is there a column vector e ∈ F n 2 , with at most p non-zero coordinates, such that He = s ? Although this problem is NP-complete [START_REF] Berlekamp | On the inherent intractability of certain coding problems (corresp.)[END_REF], in practice, the efficiency of the probabilistic algorithms devoted to solve the SD problem [START_REF] Baldi | A finite regime analysis of information set decoding algorithms[END_REF] has as a consequence that code based cryptography usually suffers from huge keys. Numerous strategies have been deployed to obtain a compact representation of the key. Among them, the use of double circulant codes [START_REF] Gaborit | Shorter keys for code-based cryptography[END_REF] leads to secure protocols with short keys.

Definition 1 An n ×n matrix is called a circulant matrix if each row is obtained from the previous one by a cyclic shift over one position to the right. In the sequel, we use some coding theory terminology. The reader may refer to [21, chapter 1] for more information on coding theory. A parity check matrix H of the (2n, n) double circulant code has the following form :

H = A M ,
where A and M are two n × n circulant matrices.

A parity check matrix of a (2n, n) double circulant code can be stored, in a compact way, using only its first row. There is no general complexity result for the SD problem where H is the parity check matrix of a random double circulant code. However, in practice, up to a small factor, the best attacks against the SD problem in this case are the same as those for random binary codes. Indeed, according to [START_REF] Gaborit | Lightweight code-based identification and signature[END_REF], when n is prime and 2 is a primitive root of Z/nZ, almost all random double circulant codes lie on the Gilbert-Varshamov bound. As a result, the SD problem is considered hard by the cryptographic community for double circulant codes.

Let y = (y 0 , . . . , y 2n-1 ) ∈ F 2n 2 a column vector and let us define y (1) (X) = y 0 + y 1 X + • • • + y n-1 X n-1 and y (2) (X) = y n + y n+1 X + • • • + y 2n-1 X 2n-1 . Given that the algebra of n × n circulant matrices over F 2 is isomorphic to the algebra of polynomials in the ring F 2 [X]/(X n -1), through the mapping ψ such that ψ(A) = a 0 + a 1 X + a 2 X 2 + • • • + a n-1 X n-1 , then the product Hy boils down to the computation of two polynomial multiplications, namely:

ψ(A) × y (1) (X) (mod X n -1)
and ψ(M ) × y (2) (X) (mod X n -1) .

BIKE [START_REF] Aragon | Bit Flipping Key Encapsulation (BIKE)[END_REF] and HQC [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF] make use of this isomorphism which maps a matrix-vector product into a polynomial multiplication in F 2 [X]/(X n -1). Due to the polynomial size (from 10000 to 60000 bits) in both protocols, it turns out that this operation has an impact on key generation, key encapsulation and key decapsulation mechanisms.

For example, in the HQC submission, key generation requires one multiplication, encapsulation requires two multiplications and decapsulation requires three multiplications. These multiplications are computed over F 2 [X]/(X N -1), with 17669 N 57637.

Moreover, the multiplications are performed with one sparse operand, while the other one is dense. Sparsedense multiplications are classically implemented using convolution approaches which is an adapted version of the schoolbook approach (for example, see Aranha et al. in [START_REF] Guimarães | Secure and efficient software implementation of qc-mdpc code-based cryptography[END_REF]).

The report [START_REF] Alagic | Status Report on the First Round of the NIST PQC Standardization Process[END_REF] mentions that side-channel resistance is a desirable security property for NIST PQC candidates. A minimum requirement for cryptographic primitives to ensure this property is to provide constant time implementations. In BIKE and HQC, some secret data are represented as sparse polynomials used as operand of a multiplication by an arbitrary polynomial, in the three steps of the protocol: key generation, encryption and decryption mechanisms. Thus, any multiplication algorithm taking advantage of the sparsity of the secret data may leak some information on it. An adversary able to exploit such source of leakage may recover information on secret data. That is why dense-dense approaches, which process the sparse operand as an arbitrary polynomial, are to be considered as mandatory.

Dense-dense multiplication over F 2 [X] has been intensively studied in the past, for different applications:

schoolbook approaches (with quadratic complexity); -Karatsuba-Offmann [START_REF] Karatsuba | Multiplication of many-digital numbers by automatic computers[END_REF] and Toom-Cook [START_REF] Bodrato | Towards optimal Toom-Cook multiplication for univariate and multivariate polynomials in characteristic 2 and 0[END_REF] subquadratic methods, with interpolation-evaluation algorithms;

-Schonhäge-Strassen [START_REF] Schönhage | Schnelle multiplikation großer zahlen[END_REF] and Fürer [START_REF] Fürer | Faster integer multiplication[END_REF] FFT based methods, and recent works (see Harvey et al. in [17,[START_REF] Harvey | Faster polynomial multiplication over finite fields using cyclotomic coefficient rings[END_REF]) showing a quasi-linear complexity in O(n log n) for integer multiplication.

One reference for the dense-dense operation is the general purpose NTL library (see [1]), which aims to provide the whole set of operations, and is based on the gf2x library [START_REF]Inria. gf2x library[END_REF] for the characteristic 2 operations. All the different approaches mentioned above are implemented. However, this library and the underlying gf2x have been designed for general purpose use and are optimized for generic operations with operand of any size.

In terms of operand size, the HQC [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF] protocol deals with polynomials whose size is fixed as a protocol parameter. The range of sizes corresponds to the one for which the Karatsuba and Toom-Cook approaches are the best in the state-of-the-art. Indeed, in the gf2x library, over F 2 [X]/(X N -1), with 17669 N 57637, the computations make use of Karatsuba or Toom-Cook (split by 3 or 4) approaches, while the threshold for FFT approaches is above 240000 bits (i.e. the degree of the polynomials). These bounds are relevant with the recent results of Harvey et al. in [START_REF] Harvey | Faster polynomial multiplication over finite fields[END_REF][START_REF] Harvey | Faster polynomial multiplication over finite fields using cyclotomic coefficient rings[END_REF].

In terms of software implementation on x86-64 platforms, until recently, the state-of-the-art was AVX2 instruction set implementations, especially using the PCL-MULQDQ instruction (see the gf2x library). This instruction performs a binary polynomial (of degree at most 63) multiplication over F 2 [X]. It returns the result stored in a 128-bit word, either an xmm 128-bit register or a same size memory storage location. The AVX2 instruction deals with 256-bit registers or memory words and performs various vectorized operation on packed operands. In 256-bit words, one can store either 32 bytes, 16 16-bit words, 8 double words or 4 quadwords (whose size is 64 bits).

In 2018, Intel announced a new instruction set extension in the so-called Icelake processor generation, which extends the AVX512 instruction set already available on some XEON processors. In particular, this architecture introduces a vectorized VPCLMULQDQ instruction, which performs up to four polynomial PCLMULQDQ multiplications, the four 128-bit results being stored in one single 512-bit word. Following this announcement, Drucker et al. [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF] proposed a software implementation (for polynomials up to degree 2 16 -1) using this instruction set, but could only experiment simulations or adapted versions of their software implementations, since no platforms and no Icelake processors were available at this time. However, they claimed up to 50% lowering of retired instructions and predicted the same drop in terms of processor clock cycle number execution. Their implementations consist in three core flows that perform schoolbook multiplication:

1. a 4 × 4 quadwords (64 bits) multiplication, written in AVX, using xmm registers and the PCLMULQDQ instruction, 2. a 4 × 4 quadwords (64 bits) multiplication, written in AVX512, using ymm registers and the VPCLMULQDQ instruction, 3. a 8 × 8 quadwords (64 bits) multiplication, written in AVX512, using zmm registers and the VPCLMULQDQ instruction.

They also mention Karatsuba multiplications for operand sizes above 256 bits. They provide the source code only for the second approach (4 × 4 quadwords multiplication, using ymm registers).

Contributions

In this work, we explore the Karatsuba and Toom-Cook multiplication construction and we identify the best combinations to be used depending on the polynomial degrees. As an illustration, we applied these results on the hqc-128 and hqc-192 multiplications of the Optimized Implementation of the HQC release, 2020/10/01 version [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF]. In this release, the multiplication implementation make use of a 2-recursive 3-split Karatsuba. We show that using a Tom-Cook-3 approach, this provides some speedup in comparison with the initial multiplication of this release. As a consequence, this has been integrated in the last official HQC release, 2021/06/06 version.

Then, in the context of AVX512 instruction set, now available since the Icelake microarchitecture processors, we propose new implementations designed for cryptographic use of polynomial multiplications over F 2 [X]. We show that the elementary multiplication construction has to be a schoolbook approach up to the 256 bit operand level, while in the state-of-the-art AVX2 context, a Karatsuba multiplication is required at the threshold of 128 bit operands.

We then implement tailor made vectorized subquadratic approaches (recursive Karatsuba and Toom-Cook-3) using the AVX512 instruction set and the vectorized VPCLMULQDQ instruction in order to improve the performances, in comparison with current state-of-theart AVX2 implementations. We compare our implementations:

with the gf2x library; -with the multiplications provided or derived from the Optimized Implementation of HQC.

Drucker et al. in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF] estimated that, by using the new AVX512 instruction set and especially the new VPC-LMULQDQ one, the retired instruction count might be divided by two, and the clock cycle number might be lowered in the same proportion. At the time they submitted their paper, there were no actual processor available implementing the instruction set extension with VPCLMULQDQ. In this paper, we checked their claims. We show that while the instruction count reduction can be overtaken, in our tests, the clock cycle number is lowered by about 39% only.

This paper is organized as follows : in Section 2 we present the Karatsuba multiplication over F 2 [X] and our AVX512 software implementation, the timing results and comparison with the implementations of Drucker et al. [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF][START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF], the gf2x library and state-of-theart AVX2 implementation; in Section 3 we present the Toom-Cook multiplication over F 2 [X] and our AVX512 software implementation, the timing results and comparison with the gf2x library and state-of-the-art AVX2 implementation; in Section 4 we present the performances obtained with the HQC protocol:

when using our AVX2 multiplications in hqc-128 and hqc-192 in the HQC release (round 3), 2020/10/01 version; -when implementing our AVX512 multiplications in the last official HQC release, 2021/06/06 version.

Finally Section 5 provides some concluding remarks. The source code of all of our implementations are available at https://github.com/arithcrypto/, in the AVX512PolynomialMultiplication repository.

Karatsuba multiplication: algorithms and implementations

The Karatsuba multiplication algorithm is the first subquadratic approach, which has been presented by Karatsuba and Offmann in [20]. This multiplication was first applied to large integers, but can be applied to polynomials. This classical approach has been extensively studied since then, and our work relies on all those previous works. Our main contribution here is the AVX512 software implementation of this approach, in order to speedup the runtime execution of multiplications over F 2 [X].

First, we review the subquadratic Karatsuba approaches for multiplication over F 2 [X]. We then present our implementations and the performance results.

Karatsuba algorithm

One wants to multiply two arbitrary polynomials of degree at most N -1, and the result is of degree at most 2 • N -2. The Karatsuba complexity applied to polynomial multiplication over F 2 [X] has been studied by Nègre and Robert in [START_REF] Nègre | Impact of Optimized Field Operations AB, AC and AB + CD in Scalar Multiplication over Binary Elliptic Curve[END_REF]. Let A and B be two binary polynomials of degree at most N -1. These polynomials are packed into an array of 64-bit words, whose size is N/64 . Let t = 2 r with r the minimum value ensuring t N/64 . Now, A and B are considered as polynomials of degree at most 64 • t -1. We reproduce the Karatsuba algorithm in Algorithm 1 Appendix A. From [START_REF] Nègre | Impact of Optimized Field Operations AB, AC and AB + CD in Scalar Multiplication over Binary Elliptic Curve[END_REF], the complexity of the recursive Karatsuba multiplication is : 8t log 2 (3) -8t XOR between 64bit words and t log 2 (3) native 64-bit multiplication. We assume that this native multiplication line 2 (denoted Mult64) is performed using a single processor instruction: this is the case of the Intel Cores i3, i5 and i7 and above.

There are variants of these approaches, splitting the operands in all number of parts, and using an elementary multiplication which can be all sort of Karatsuba multiplication for example. These variants have been extensively studied in Weimerskirch and Paar in [START_REF] Weimerskirch | Generalizations of the karatsuba algorithm for efficient implementations[END_REF]. Algorithms 2 and 3 present the 3-way and 5-way split Karatsuba (see Appendix A).

In Table 1, we remind the complexity of recursive Karatsuba multiplication (t is the size of the operands in 64-bit words).

Apart the Schoolbook, which presents the worst complexity in terms of elementary multiplications, one can verify that the Karatsuba approaches are ordered in growing complexity for equivalent sizes.

AVX512 Implementation

We propose here a little survey of the possible approaches for AVX512 implementations. Our goal is to review the state-of-the-art (to our knowledge) and possibly propose improvements. We started to evaluate the implementation approaches from [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF], which were based on the schoolbook algorithm for the 256 bit and 512 bit operand sizes. Their main claim is a reduction of 50% of the instruction count for the kernels they presented in this paper, and while an actual processor were not available at the time of their work, they evaluated a similar improvement in terms of clock cycle number for the computation time. In a more recent work, Drucker et al. in [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] proposed a new approach based on the Karatsuba algorithm. In our evaluation, this last work outperforms the first of [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF]. In this section, we briefly review the main feature of the AVX512 instruction set, and then present the most interesting approaches for elementary multiplication for the 512 bit operand size, from [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF]

8 × t log 2 (3) Karat3 6 × ((t/3) log 2 (3) ) 48 × ((t/3) log 2 (3) ) Karat5 15 × ((t/5) log 2 (3) ) 120 × ((t/5) log 2 (3) )
Table 1: Complexity of the Karatsuba's multiplication variants

AVX512 instruction set and special features

The AVX512 are 512 bit extensions to the SIMD (Single Instruction Multiple Data) 256 bit AVX (Advanced Vector Extension) instructions for x86 instruction set architecture. This was proposed by Intel since 2013 and consists of multiple extensions. In the AVX512 processors, in addition to the general purpose 64 bit registers, larger registers are also available in order to perform vectorized instructions. These registers are of type xmm : of size 128 bits, i.e. containing two quadwords, thus denoted {a 1 , a 0 }, a 0 and a 1 being the quadwords in register a; ymm : of size 256 bits, i.e. containing four quadwords, thus denoted {a 3 , a 2 , a 1 , a 0 }, a 0 , a 1 , a 2 and a 3 being the quadwords in register a; zmm : of size 512 bits, i.e. containing eight quadwords, thus denoted {a 7 , a 6 , a 5 , a 4 , a 3 , a 2 , a 1 , a 0 }, a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 and a 7 being the quadwords in register a;

The number of registers of each type is 32.

The elementary polynomial multiplication:

Let us first remind how the PCLMULQDQ instruction works.

The intrinsic available by including the immintrin.h file is :

__m128i _mm_clmulepi64_si128 (__m128i a, __m128i b, const int imm8)
The quadwords a i and b i represent binary polynomials of degree at most 63. The PCLMULQDQ instruction returns the results in an xmm register, that is a j × b i , of degree at most 126. The selection of i and j, i.e. the corresponding quadword of the operand is made according to the value of imm8 :

1. imm8 = 0x00 : i = 0, j = 0 → PCLMULQDQ returns a 0 × b 0 ; 2. imm8 = 0x01 : i = 0, j = 1 → PCLMULQDQ returns a 1 × b 0 ; 3. imm8 = 0x10 : i = 1, j = 0 → PCLMULQDQ returns a 0 × b 1 ; 4. imm8 = 0x11 : i = 1, j = 1 → PCLMULQDQ returns a 1 × b 1 ;
The VPCLMULQDQ instruction now available on Icelake and above platforms has the following intrinsic:

__m512i _mm512_clmulepi64_epi128 (__m512i a, __m512i b, const int Imm8)
This instruction computes in parallel 4 pclmul multiplications, i.e. carryless multiplications of binary polynomials of degree at most 63, stored in 4 quadwords in 512 bit registers, as seen above. Thus, four 128 bit results are stored in the zmm register as follows:

{a 6+j × b 6+i 128 bits , a 4+j × b 4+i 128 bits , a 2+j × b 2+i 128 bits , a j × b i 128 bits

}

The selection of i and j is made as above according to the value of Imm8.

We now examine the implementation of four multiplications using this instruction set and these registers: the mul512 version, from Drucker et al. in [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF], our new karat_1_512 using the mul128x4 procedure from [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] and our new karat_mult_1_512_SB using a schoolbook mul128x4 procedure, and the full schoolbook 512 bit implementation, as suggested in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF]. We provide the detailed source code of the karat_mult_1_512_SB, and corresponding explanations in Appendix C, while the source code of the other approaches are available in the github repository, as mentioned in the Introduction.

We chose not to present here the schoolbook approaches of [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF], because in our tests, these versions are outperformed by the others. Likewise, we also implemented 256 bit kernels using schoolbook (SB256) and Karatsuba at the 256 bit level Karat256 with the AVX512 and VPCLMULQDQ instruction. These versions are also outperformed by the others, in particular by the 8x8 multiplication of [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF].

Nevertheless, we give an overview on these multiplications and their performances Appendix B and Appendix F, and the source code is also provided in the github repository.

mul512 version, from Drucker et al. in [11]

In [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF], Drucker et al. present a multiplication of 1024 bit operands. This multiplication is computed as follows:

-The mul1024 is a Karatsuba wrapper which calls three mul512 multiplications, along with a classical Karatsuba reconstruction. This implementation is similar to the AVX2 equivalent implementations except the register size. -The mul512 is a Karatsuba multiplication which splits in four parts the 512 bit operands. They use a four 512 bit word table, storing the 5 elementary xored operands in addition to the operands themselves, for a total of 9 pairs of 128 bit operands. After this step of operand preparation, the mul512 procedure calls three times a mul128x4 function in order to compute the 9 elementary 128 bit operand multiplications. -Finally, the mul128x4 procedure performs four 128 bit operand multiplications in parallel, using the VPCLMULQDQ instruction. This instruction is called three times, corresponding to the Karatsuba construction using 128 bit operands, split in two 64 bit words.

One may notice that the mul512 procedure invokes 9 times the VPCLMULQDQ instruction, that is 36 elementary 64 bit operand multiplications, while using only 3 3 = 27 out of them due to the Karatsuba multiplication. We refer the reader to [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF]Appendix B], for a complete and detailed explanation of the source code.

2.2.4

Our new karat_mult_1_512 using the mul128x4 procedure, from [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] Starting from the previous implementation of [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF], our goal is to check the difference between the Karatsuba and the schoolbook approach at the 256 bit operand level. For this sake, we modified the mul512 procedure into classic Karatsuba construction, which split in two parts the 512 bit operands. Now, three elementary 256 bit schoolbook multiplications are computed. These multiplications invoke one single call to the mul128x4 procedure from [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF]. In our version, the code of the mul128x4 procedure has been manually inlined.

The total number of VPCLMULQDQ instructions in this 512 bit multiplication is 9, the same as previously. However, this procedure now makes use of all 36 elementary 64 bit multiplication computed, while slightly simplifying the final reconstruction at the 256 bit level. This version is called karat_mult_1_512, and its source code can be found in the github repository of the paper.

Our new karat_mult_1_512_SB using a schoolbook mul128x4 procedure

This configuration implements the schoolbook algorithm at the 128 bit and 256 bit multiplication levels. Our goal now is to check which algorithm between Karatsuba and schoolbook is the best at the 128 bit operand level. Indeed, the instruction count is lower, while making use of one more VPCLMULQDQ instruction. The latency and throughput of this instruction is higher than conventional instructions. However, the vectorized version changes this by performing simultaneously four elementary 64 bit operand multiplications.

This approach now uses 3 × 4 2 = 48 elementary 64 bit operand multiplications in total. This version is called karat_mult_1_512_SB and is presented in details in Appendix C page 17.

Full schoolbook 512 bit multiplication SB512

From the description of Drucker et al. in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF], we also wrote a full schoolbook approach at the 512 bit size level.

This version now uses 4 × 4 2 = 64 elementary 64 bit operand multiplications in total, and is called SB512.

Instruction counts for the four 512 bit multiplication versions

The instruction count for both configurations of mul128x4 is shown Table 2.

Instruction counts for mul512, SB512, karat_mult_1_512 and karat_mult_1_512_SB are shown Table 3. For the Drucker et al. [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] version, we count three times the mul128x4 instruction count plus the mul512 instructions.

These instruction counts give an overview of the potential differences between the schoolbook and Karatsuba approaches at different levels. The threshold for the use of Karatsuba algorithm is at the lowest level in AVX2 implementation. However, the vectorized VPCLMULQDQ instruction performs 4 multiplications at a time, with similar latency and throughput. Thus, in the context of AVX512, the good level for this threshold is 256 bit multiplication.

As expected (see Table 3), the instruction count is the lowest for the karat_mult_1_512_SB, corresponding to the schoolbook approach until the 256 bit operand size. The mul512 of [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] instruction count is the greatest among the three approaches, due to the more complex final reconstruction of each Karatsuba recursion step. We present Table 4 the results for the first three approaches described above, and compare them to the classic AVX2 implementations from the gf2x library [START_REF]Inria. gf2x library[END_REF] and the one from the HQC submission [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF]. The main conclusions are as follows:

among the AVX512 implementations, the best version is the one based on the karat_mult_1_512_SB, i.e. a schoolbook algorithm applied to the 128 bit and 256 bit levels, and the Karatsuba approach at the highest level (512 bit operands); -in comparison with the AVX2 implementation, this version achieves a retired instruction number divided by roughly three, and speedups (in terms of clock cycle number) are from 29.5% to up to 39.1% (KaratRec, size 131072 bit), due to the lower IPC of the AVX512 instruction set.

The gf2x library performances are slightly lower than the ones of the other implementations. The results shown Table 4 allows to evaluate the Instructions per cycles (IPC), which is the ratio between the retired instructions and the clock cycle number. The graph in Figure 1 shows that while the IPC of the gf2x library multiplications is about 4.5, the IPC of the AVX2 multiplications is nearly 3.0 and the IPC of our AVX512 multiplications (with karat_mult_1_512_SB) is about 1.5. In the gf2x case, the software implementation makes use of AVX2 instruction set. However, this all-purpose library includes some wrappers and tests, and especially offers the possibility to tune each operand size. This is costly, but is written with conventional instructions. This is the most likely explanation of the high IPC, while the clock cycle number is a little worse that the one of the AVX2 software implementation.

We compare also the IPC of our AVX512 multiplication (with karat_mult_1_512_SB) and the one using the mul512 of Drucker et al. from [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF]. In this last case, the IPC is greater than the one obtained with our implementation. While the clock cycle numbers of our implementations are about 9% lower than the ones using the mul512, The conclusion of this comparison is that the AVX512 multiplication implementation presents an elementary 64 bit multiplication cost relatively low in comparison with the AVX2 situation. In this case, the non vectorized multiplication cost leads to a Karatsuba application at the 128 bit level. Insofar as the vectorized instructions equivalently divides by four the latency and the throughput of the elementary multiplication at the 64 bit level, it becomes more interesting to apply the schoolbook approach at the 128 bit level and also at the 256 bit level. In our tests, we saw that at the 512 bit level, the Karatsuba approach becomes again the best one.

On the hardware point of view, our tests show another aspect of this AVX512 implementation case, which differs from the initial estimations given in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF], concerning the potential speedups brought by the AVX512 instruction set and VPCLMULQDQ instruction. Indeed, these results make clear that, in terms of processor hardware, at the microarchitectural level, the AVX512 instruction set has not yet been implemented with the same integration level as the one of the other instruction sets (conventional and AVX2), at least on our platform. The Intel documentation does not provide a lot of details [START_REF]Intel. Intel® 64 and ia-32 architectures software developer manuals[END_REF], however, one can assume that the hardware features are not homogeneous with the AVX2 equivalents. This also means that future Intel processor generations might improve the IPC of the AVX512 instruction set and potentially decrease the clock cycle number of our implementation. Indeed, between our AVX2 and AVX512 multiplication implementation, the retired instruction count is divided by nearly three. If the processor manufacturer improves the IPC of the AVX512 and AVX2 instruction sets, to get closer to the one of the conventional instruction set IPC, this means that one may observe speedups with our software implementations on future platforms.

Karat3 and Karat5 implementations

These multiplications (Algorithms 2 and 3, Appendix A page 15) make use of the previous KaratRec multiplication (Algorithm 1) as elementary multiplications. The vectorized versions are again implemented using AVX2 and AVX512 instruction set. This leads to different sizes: while the Karat3 multiplication has an operand size which is three times the one of the elemen- Fig. 1: Instructions per cycle comparison for Algorithm 1 tary multiplication (for example 512,1024, 2048... bits), the Karat5 has an operand size which is five times the one of the elementary recursive Karatsuba multiplication. Thus, depending on the context, one may choose the most appropriate version according to the operand sizes. We present the performance results (# clock cycles) in Table 5.

The Karat3 or Karat5 multiplications can also use themselves as elementary multiplications. This leads to four extra combinations for multiplication. We present the performance results (# clock cycles) in Table 6. One may notice that the Karat3(Karat5), i.e. a Karat3 multiplication using a Karat5 elementary multiplication, deals with the same operand sizes that the Karat5(Karat3). However, we give the results for both versions, and observe that they present very close performances.

In all this experimentation, we only used, as the most elementary multiplication, the karat_mult_1_512_SB version, presented Section 2.2.5 page 6. The platform and experimentation process are the same as previously (see Appendix E page 19).

Conclusion

In this section, we presented our AVX512 implementation of recursive Karatsuba multiplication over F 2 [X] for polynomial of degree at most 131071.

Our implementations show that the best AVX512 approach is a 512 bit kernel, using a schoolbook algorithm at the 128 and 256 bit level, and Karatsuba at the highest level, that is for the 512 bit size operands.

Used in AVX512 recursive Karatsuba multiplications of greater sizes, and in comparison with the AVX2 software implementation of [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF] a retired instruction number divided by roughly three, and speedups (in terms of clock cycle number) are from 29.5 % to up to 39.1 % (KaratRec, size 131072 bit), due to the lower IPC of the AVX512 instruction set.

The same achievements have been reached in all the Karatsuba variants (split in 3, 5 parts, and combinations).

Toom-Cook multiplication over F 2 [X]

In this section, we present the implementation issues of the Toom-Cook 3-5 approach applied to multiplication over F 2 [X], especially with the AVX512 instruction set. We refer the reader to Appendix D page 18 for a general presentation of the Toom-Cook approach.

Toom-Cook multiplication complexity

There are several way to split the operands in the Toom-Cook approach:

-The Toom-Cook 3-5, which splits the operands in 3 parts, and involves 5 elementary multiplications; -The Toom-Cook 4-7, which splits the operands in 4 parts, and involves 7 elementary multiplications; -The Toom-Cook 5-9, which splits the operands in 5 parts, and involves 9 elementary multiplications;

We present Table 7 the general results versus the polynomial degrees, and the corresponding operand size in 64 bit words. As log 2 (3) > log 3 (5) > log 4 (7) > log 5 [START_REF] Brent | Faster multiplication in gf(2)[x[END_REF], this implies that the multiplication number is decreasing while increasing the split. However, the hidden constant in the O notation is increasing, due to the more complex interpolation phase. This leads to a "gray zone" in which the different algorithms are very close to each other and this gray zone resides in our range for code based cryptographic protocols. In the gf2x implementation, a specific selection of the algorithm is made depending on each of the operand size (which can be different in this general purpose library).

Toom-Cook multiplication, implementation issues

Let us recall that in order to multiply two binary polynomials A and B of degree at most N -1, we consider them as polynomials of degree at most 64t -1 where t = 3n and n ensures t N/64 . We now present how to choose n.

In the evaluation phase, the elementary products do not have the same operand size: C(0), C(1) and C(∞) have operands of degree at most 64n-1, while C(x) and C(x + 1) have operands of degree at most 64n + 2 • w -1 (see Appendix D page 18). To take into account this characteristic, if the size of the elementary product is known, one has to set operands of size n + 2w/64 64bit words, padding with zeros in order to use the same elementary product. Now, we can specify the value of n mentioned Appendix D: n must be the minimum value ensuring 3n

N/64 such that n + 2w/64 is the size of the elementary multiplications computed during the evaluation phase. Thus, n + 2w/64 is either a power of 2 in case of a Karatsuba multiplication (as seen section 2.1), or a value which complies with another Toom-Cook multiplication.

In the interpolation phase, since the divisions by x and (x+1) are exact, they can be implemented using the trick presented by Quercia and Zimmermann in [START_REF] Zimmermann | Irred-ntl patch[END_REF] and [START_REF] Quercia | Irred-ntl source code[END_REF]. One takes advantage of the size of the polynomial to replace these divisions by a one word right shift for the division by x, and by a special multiplication by (x+ 1) -1 mod X d , d > degree of C(x) and d ≡ 0 mod w.

In the second case, to evaluate (x + 1) -1 mod X d = (X w + 1) -1 mod X d , notice that:

(X w + 1) -1 mod X n = d/w-1 i=0 X w•i . Example 1.
Here is a small toy example with polynomial over F 2 [X] : let us divide by X + 1 polynomial whose degree is less than 8. One has (X + 1) -1 mod X 8 = X 7 + X 6 + X 5 + X 4 + X 3 + X 2 + X + 1. Now, ∀P (X) such as (X + 1)|P (X) and degree(P )< 8, Q(X) = P (X)/(X +1) is computed as P (X)•(X +1) -1 mod X 8 . Since this division is exact, the result is the exact quotient. Let us set P (X) = X 7 + X 5 + X 4 + X :

Q(X) = (X 7 + X 5 + X 4 + X) • (X + 1) -1 mod X 8 = (X 7 + X 5 + X 4 + X) •(X 7 + X 6 + X 5 + X 4 + X 3 + X 2 + X + 1) mod X 8 = X 6 + X 5 + X 3 + X 2 + X.
The vectorized implementation, while using 256-bit instructions, uses both values w = 64 or w = 256. We can even use w = 512 in the case of AVX512 platforms. Notice that the division by x + 1 is cheaper in case of w = 256 or 512, while n is slightly greater with w = 64.

To illustrate some Toom-Cook use cases, let us present some examples.

Example 2. In order to use elementary Karatsuba multiplications, and with w = 64, let us consider n such that n+2w/64 = 2 8 . Thus, one has 64•n+2•w = 8192, and this elementary multiplication proceeds polynomials of degree at most 8191. We then have n = 254, t = 3 • n = 762, thus building a Toom-Cook multiplication which can multiply polynomials of degree at most the first split computes polynomials of degree at most 23679; we build the operands for the evaluation phase, whose size is up to 24192; we then call the Toom-Cook procedure of the previous example, which makes use of Karatsuba multiplications of size 8192.

Example 4. We now build a Toom-Cook multiplication, setting w = 512 using a Karat5 elementary multiplication of size 20480. One now chooses n such that 64n+2•512 = 20480. This leads to n = 304 and t = 912. This multiplication proceeds polynomials of degree at most 58368. This multiplication fits with the hqc-256 protocol (see [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF]), whose parameter N = 57669, and for an AVX512 implementation. This can be adjusted for other sizes. To deal only with word shifts, this implies on the operand size that:

elementary recursive Karatsuba multiplications are based on PCLMULQDQ or VPCLMULQDQ 64-bit elementary multiplications, whose size corresponds to the ones previously seen (KaratRec, Karat3, or Karat5...); -the computation of C(x) and C(x + 1) needs construction of operands whose size has to be the elementary multiplication size; -the split has to take into account the size of the elementary multiplication operand, by diminishing the size of the split by two words.

Toom3Mult implementations

Three Toom3Mult versions have been implemented (see The speedup of our AVX512 implementation in comparison with the AVX2 one is as follows:

-Toom3 based on KaratRec (Toom3Mult(KaratRec)):

the speedup starts from 33.5% for the 48768 bit operand size up to 36.7% for 97920 bit operand size; -Toom3 based on Karat3 (Toom3Mult(Karat3)): the speedup starts from 30.4% for 36480 bit operand size up to 33.3% for the 73344 bit operand size; -Toom3 based on Karat5 (Toom3Mult(Karat5)): the speedup starts from 32.4% for the 14976 bit operand size up to 33.9% for 122496 bit operand size.

Toom3 vs Karat3 comparison

The Toom-Cook multiplication presented above (Toom3) and the Karat3 multiplication (see Algorithm 2 page 15) look similar, both splitting the operands in three parts. However, the Toom3 needs 5 elementary multiplications while the Karat3 requires 6. In the gf2x library [START_REF]Inria. gf2x library[END_REF], the Toom-3 multiplication is used above the threshold of 21 64-bit words, i.e. polynomial of degree 1343.

In our AVX2 and AVX512 implementations, we compare the clock cycle numbers of multiplications using the same elementary Karatsuba multiplication. Consequently, the operand size in the Toom3 case is slightly lower. In Table 9, the clock cycle numbers of the Toom3 multiplication are lower for the considered sizes, and the speedup starts from 8% for the smallest (about 6000 bit operand size) up to 11.7% for the bigger sizes, the maximum potential speedup being theoretically 16.7%. Indeed, the costliest interpolation and reconstruction phase of the Toom3 approach lowers the speedup for the smallest sizes.

Conclusion

In this section, we presented our AVX512 implementation of Toom-Cook multiplication over F 2 [X] for polynomial of degree at most 122496, based on various Karatsuba version elementary AVX512 multiplications (KaratRec, Karat3 and Karat5).

The speedup between AVX2 and our AVX512 implementation is again up to nearly 37% (97920 bit size, Toom3Mult(KaratRec)), as it has already been observed with our AVX512 Karatsuba multiplications. Thus, the same remark can be done about the potential in terms of future results.

HQC multiplications

In this section, we present the application of our AVX2 and AVX512 Toom-Cook multiplications in the context of the HQC protocol.

The sizes of the HQC multiplications are (see [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF]):

hqc-128 : PARAM_N = 17669; hqc-192 : PARAM_N = 35851; hqc-256 : PARAM_N = 57637.

In the NIST submission [2, updated submission package (round 3), 2020/10/01], the hqc-128 and the hqc-192 implementations make use of the Karat3 multiplication based on elementary Karat3 multiplications, while the hqc-256 is a Toom-Cook multiplication (3 part operand split) based on a Karat5 elementary multiplication. The reader may notice that the operand size is different in the hqc-256 case, in comparison with the Toom3Mult(Karat5) version presented in Table 8 page 12, whose operand size is 61056 bits, while the AVX2 operand size is 59904 bits, and the AVX512 one is 58368 bits, see Table 10. This is due to the word size considered in the Toom-Cook implementation: the first version of table 8 has a word size w = 64 bits, while the AVX2 version has w = 256 and we chose w = 512 for the AVX512 version. The new sizes fit the hqc-256 PARAM_N = 57637, allowing a slightly better performance.

The platform and test procedure are the ones described Appendix E page 19. We kept the compiler flags used in the NIST submission [2, updated submission package (round 3), 2020/10/01]: -O3 -funroll-all-loops -flto -march=tigerlake for AVX512 versions; -O3 -funroll-all-loops -flto -mavx -mavx2 -mbmi -mpclmul for AVX2 versions.

In Table 10, we sum-up the performances of the respective multiplications, from Tables 8 page 12 and 9 page 12, along with the specific hqc-256 multiplications results, in the context of HQC protocol.

HQC implementation performances

This leads to the performances and also the speedups brought by our work in comparison with the NIST release [2, updated submission package (round 3), 2020/10/01], provided in Table 11.

As a consequence of this work, the AVX2 Toom3 multiplication based on Karat3 is now part of the last release of HQC (2021/06/06).

Table 12 provides the HQC performances when using our AVX512 multiplications in the 2021/06/06 release. One may notice that column one of this Table is different of the corresponding one of Table 11 because other parts of the HQC source code has been updated between both releases.

A maximum 11.8% speedup for the HQC implementation is reached with the hqc-256 Keygen.

Conclusion

In this paper, we considered the software AVX512 implementation of polynomial multiplication over F 2 [X], using the vectorized 64-bit polynomial multiplication instruction VPCLMULQDQ. We studied the different combinations of schoolbook/Karatsuba constructions for the kernels up to 512 bit operands. We then implemented two different approaches: one based on the Karatsuba subquadratic approach and the other on the Toom-Cook approach. These implementations are competitive in comparison with state-of-the art general purpose library, HQC submissions, and other AVX512 software implementations of [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF][START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF]. While the retired instruction count is divided by roughly three compared to the corresponding AVX2 implementations, we achieved a speedup up to nearly 40%, in terms of clock cycle numbers.

We implemented our approaches in the HQC protocol by patching the NIST submission released in october 2020, in order to experiment the potential benefits, and this leads to speedups up to 11.8% (hqc-256 Keygen.).

All the implementations of this work are available on github 1 .

Funding This work has been partially funded by TPM Metropol (AAP2020-IPOCRAS project).

Algorithm 3 Karat5(A,B,t), from [START_REF] Weimerskirch | Generalizations of the karatsuba algorithm for efficient implementations[END_REF] Require: A and B on t = 5 × 2 r computer words. Ensure: R = A × B 1: // Split in five parts of word size t/5. 2: A = A 0 + x 64t/5 A 1 + x 2×64t/5 A 2 + x 3×64t/5 A 3 + x 4×64t/5 A 4 3: B = B 0 + x 64t/5 B 1 + x 2×64t/5 B 2 + x 3×64t/5 B 3 +

x 4×64t/5 B 4 4: // Recursive multiplications 5: R 0 ← KaratRec(A 0 , B 0 , t/5) These lines define the constant indexes for the _mm512_permutexvar_epi64 and _mm512_permutex2var_epi64 instructions. These instructions are explained on the fly. The _mm512_broadcast_i64x4(*A256) instruction duplicates the 256 bits of *A256 in the A512 register, the same for the *B256.

+ (R 0 + R 1 + R 01 )X 64t/5 + (R 0 + R 1 + R 2 + R 02 )X 2×64t/5 + (R 0 + R 1 + R 2 + R 3 + R 03 + R 12 )X 3×64t/5 + (R 0 + R 1 + R 2 + R 3 + R 4 + R 04 + R 13 )X 4×64t/5 + (R 1 + R 2 + R 3 + R 4 + R 14 + R 23 )X 64t + (R 3 + R 2 + R 4 + R 24 )X 6×64t/5 + (R 3 + R 4 + R 34 )X 7×64t/5 + R 4 X 8×64t/
The _mm512_permutexvar_epi64 (idx_b, tmp) spreads the 64 bit words following the index idx_b. This allows to shuffle the 64 bit words of *B256 in the B512 register, in order to prepare the elementary multiplications.

We thus have:

A512 ← {a 3 , a 2 , a 1 , a 0 , a 3 , a 2 , a 1 , a 0 } tmp ← {b 3 , b 2 , b 1 , b 0 , b 3 , b 2 , b 1 , b 0 } B512 ← {b 1 , b 0 , b 3 , b 2 , b 3 , b 2 , b 1 , b 0 } R0_512=_mm512_clmulepi64_epi128 ( A512 , B512 , 0 x00 ) ; R1_512=_mm512_clmulepi64_epi128 ( A512 , B512 , 0 x10 ) ; R2_512=_mm512_clmulepi64_epi128 ( A512 , B512 , 0 x01 ) ; R3_512=_mm512_clmulepi64_epi128 ( A512 , B512 , 0 x11 ) ;
We now compute all the elementary 64 bit multiplications, providing all the 128 bit results as follows:

R0 ← {a 2 × b 0 , a 0 × b 2 , a 2 × b 2 , a 0 × b 0 } R1 ← {a 2 × b 1 , a 0 × b 3 , a 2 × b 3 , a 0 × b 1 } R2 ← {a 3 × b 0 , a 1 × b 2 , a 3 × b 2 , a 1 × b 0 } R3 ← {a 3 × b 1 128bits , a 1 × b 3 128bits , a 3 × b 3 128bits , a 1 × b 1 128bits } tmp = _mm512_permutex2var_epi64 ( R0_512 , idx_1 , R3_512 ) ;
The tmp register now contains all the a i × b i elementary products coming form the R0_512 and R3_512 registers:

tmp ← {a 3 × b 3 , a 2 × b 2 , a 1 × b 1 , a 0 × b 0 }
It remains now to compute the middle part of the result to be added to tmp, in order to get the final result. The middle register now contains the addition (XOR) between R1_512 and R2_512, reordered with the _mm512_permutexvar_epi64:

middle

← {a 0 b 3 ⊕ a 3 b 0 , a 2 b 3 ⊕ a 3 b 2 , a 1 b 2 ⊕ a 2 b 1 , a 1 b 0 ⊕ a 0 b 1 } tmp ^= _mm512_permutex2var_epi64 ( middle , idx_4 , idx_b ) ; tmp ^= _mm512_permutex2var_epi64 ( middle , idx_5 , idx_b ) ;
The tmp register is added (XOR) with the elementary products of the middle register, and nearly contains the result, except some of the products of the middle part: The remaining products of the middle part to be added with tmp are put in place in the middle register:

tmp ← a3 × b3 a2 × b2 a1 × b1 a0 × b0 ⊕ a2b3 ⊕ a3b2 a1b2 ⊕ a2b1 a1b0 ⊕ a0b1 ⊕ a0b3 ⊕ a3b0 m i d d l e = _mm512_permutex2var_epi64 ( R0_512 , idx_6 ,
middle ← {0x0U L 128bits , a 1 b 3 ⊕a 3 b 1 , a 0 b 2 ⊕a 2 b 0 , 0x0U L 128bits } * Out = tmp^m i d d l e ; }
Out gets the final reconstruction :

Out ← a3 × b3 a2 × b2 a1 × b1 a0 × b0 ⊕ a2b3 ⊕ a3b2 a1b2 ⊕ a2b1 a1b0 ⊕ a0b1 ⊕ a1b3 ⊕ a3b1 a0b2 ⊕ a2b0 ⊕ a0b3 ⊕ a3b0
// t h i r d m u l t i p l i c a t i o n 256 : SASB . . . cm ^= m i d d l e ^c l ^ch ;

The register ch now contains the 512 bit result of ah×bh.

The result cm is directly added (XOR) to the other results cl and ch in order to prepare the final Karatsuba reconstruction, and cm now contains the 512 bit result of (sa×sa)⊕cl⊕ch: This ends the computation, the final lines stores the result: the 512 least significant bits in the memory place C[0], and the most significant bits in C [1].

Appendix D Toom-Cook multiplication general algorithm

Several approaches to multiply two arbitrary polynomials over F 2 [X] of degree at most N -1, using the Toom-Cook algorithm, have been presented by Bodrato in [START_REF] Bodrato | Towards optimal Toom-Cook multiplication for univariate and multivariate polynomials in characteristic 2 and 0[END_REF], Brent et al. in [START_REF] Brent | Faster multiplication in gf(2)[x[END_REF], and software implementations have been provided by Quercia and Zimmermann, in the context of the ntl and the gf2x library, see [START_REF] Zimmermann | Irred-ntl patch[END_REF] and [START_REF] Quercia | Irred-ntl source code[END_REF]. Let A and B be two binary polynomials of degree at most N -1. These polynomials are packed into an array of 64-bit words, whose size is N/64 . Let t = 3n with n a value ensuring t N/64 . Now, A and B are considered as polynomials of degree at most 64 • t -1. We discuss the value of n in section 3.2.

A and B are split in three parts. One wants now to evaluate the result

C = A • B with A = a 0 + a 1 • X 64n + a 2 • X 2•64n ∈ F 2 [X], B = b 0 + b 1 • X 64n + b 2 • X 2•64n ∈ F 2 [X],
(of maximum degree 64t -1, and a i , b i of maximum degree 64n -1) and,

C = c 0 + c 1 • X 64n + c 2 • X 2•64n + c 3 • X 3•64n + c 4 • X 4•64n of maximum degree 6 • 64n -2.
The "word-aligned" version evaluates the polynomial for the values 0, 1, x = X w , x + 1 = X w + 1, ∞, w being the word size, typically 64 in modern processors. Furthermore, on Intel processors, one can set w = 256 to take advantage of the vectorized instruction set AVX-AVX2, and even w = 512 (AVX512 extension), at the cost of a slight operand size reduction. For the evaluation phase, one has:

C(0) = a 0 • b 0 C(1) = (a 0 + a 1 + a 2 ) • (b 0 + b 1 + b 2 ) C(x) = (a 0 + a 1 • x + a 2 • x 2 ) • (b 0 + b 1 • x + b 2 • x 2 ) C(x + 1) = (a 0 + a 1 • (x + 1) + a 2 • (x 2 + 1))• (b 0 + b 1 • (x + 1) + b 2 • (x 2 + 1)) C(∞) = a 2 • b 2
The implementation of this phase is straightforward, providing that the multiplication a i •b i is either another Toom-Cook or Karatsuba multiplication. Notice that the multiplications by x or x 2 are virtually free word shifts.

For the interpolation phase, one has the following equations:

C(0) = c 0 C(1) = c 0 + c 1 + c 2 + c 3 + c 4 C(x) = c 0 + c 1 • x + c 2 • x 2 + c 3 • x 3 + c 4 • x 4 C(x + 1) = c 0 + c 1 • (x + 1) + c 2 • (x 2 + 1) +c 3 • (x 3 + x 2 + x + 1) + c 4 • (x 4 + 1) C(∞) = c 4
The matrix associated to this system of equations is given by:

M =       1 0 0 0 0 1 1 1 1 1 1 x x 2 x 3 x 4 1 x + 1 x 2 + 1 x 3 + x 2 + x + 1 x 4 + 1 0 0 0 0 1      
and one has :

M -1 =        1 0 0 0 0 (x 2 +x+1) (x 2 +x) 1 1/x 1 x+1 x 2 + x 0 1 x 2 +x 1 x+1 1/x x 2 + x + 1 1 x 2 +x 1 x 2 +x 1 x 2 +x 1 x 2 +x 0 0 0 0 0 1       
Finally, the interpolation phase gives : The compiler is gcc version 10.2.0, the compiler options are as follows:

c 0 = C(0) c 1 = (x 2 + x +
-O3 -g -march=tigerlake -funroll-all-loops -lm -lgf2x.

We kept the -funroll-all-loops option though it does not provide significant improvements. We follow the same kind of test procedure that the one described in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF] :

the Turbo-Boost® is deactivated during the tests; -1000 runs are executed in order to "heat" the cache memory; -one generates 50 random data sets, and for each data set the minimum of the execution clock cycle numbers over a batch of 1000 runs is recorded; -the performance is the average of all these minimums; -this procedure is run on console mode, to avoid system perturbations, and obtain the most accurate cycle counts.

The clock cycle counter is rdtsc and the instruction counter is rdpmc with the corresponding selection. The results for the smallest sizes (i.e. 256 bit and 512 bit operand sizes) are not very reliable since rdtsc and rdpmc are not serializing instructions (see [START_REF]Intel. Intel® 64 and ia-32 architectures software developer manuals[END_REF]). For such sort of small functions, we wanted to avoid the insertion of a costly serializing instruction as cpuid, while the instruction count and the clock cycle number may be less than 20. We chose not to present them. The first size considered is 1024 bits, i.e. binary polynomial of degree at most 1023 operands.

Appendix F Instruction count and performances

F.1 Instruction count comparison

In Table 13, we provide the comparison between the instruction count of our schoolbook and Karatsuba versions. Moreover, we compare this two approaches with the current state-of-the-art AVX2 reference. Such an AVX2 implementation can be found in the source code of the optimized version of HQC [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF]. It uses the AVX2 instruction set and the non vectorized PCLMULQDQ instruction. Finally, we also put in Table 13 the instruction number of the assembly source code for the same multiplication presented by Drucker et al. in [START_REF] Drucker | Fast multiplication of binary polynomials with the forthcoming vectorized vpclmulqdq instruction[END_REF]. Here are some comments on these results:

-The best version is our implementation of the schoolbook approach, dividing by more than 2 the instruction number in comparison with the state-of-the-art AVX2 implementation. -Our Karatsuba approach presents more instructions but only 3 VPCLMULQDQ instead of 4 for the schoolbook version. Thus, the performance comparison may vary according to the latency and throughput of the instructions. -Drucker et al.'s version has 8 VPCLMULQDQ instructions and a larger instruction number (31, instead of 19 for our implementation of the schoolbook approach). This is due to the fact that they only use 2 elementary 64 bit multiplications per VPCLMULQDQ instruction (ymmm version of the instruction), while we use 4. This also implies more XOR's in their case.

F.2 Performances for the 256 bit level kernels

We present Table 14 the performances of the AVX512 Karatsuba multiplications using the 256 bit kernels presented above. We also include the results of the multiplications using our 8x8 SB-512 kernel.

a 1

 1 1 . . . a n-2 a n-1 a n-1 a 0 . . . a n-3 a n-2 a 2 . . . a n-1 a 0

Definition 2 2 .

 22 Let k, n ∈ N, an (n, k) linear code C over F 2 is a k dimensional subspace of F n Definition 3 A parity check matrix of an (n, k) linear code C is an (n -k) × n matrix H over F 2 such that H t c = 0 iff c ∈ C.Definition 4 A (2n, n) double circulant code C is a linear code such that : (c 0 , c 1 , . . . , c n-2 , c n-1 , c n , c n+1 , . . . , c 2n-2 , c 2n-1 ) ∈ C ⇓ (c n-1 , c 0 , . . . , c n-3 , c n-2 , c 2n-1 , c n , . . . , c 2n-3 , c 2n-2 ) ∈ C .

5 22 :

 22 return (R) __m512i idx_1=(__m512i ) { 0 x0UL , 0 x1UL , 0 x8UL , 0 x9UL , 0x2UL , 0 x3UL , 0 xaUL , 0 xbUL } ; __m512i idx_2=(__m512i ) { 0 x0UL , 0 x1UL , 0 x6UL , 0 x7UL , 0x2UL , 0 x3UL , 0 x4UL , 0 x5UL } ; __m512i idx_3=(__m512i ) { 0 x0UL , 0 x1UL , 0 x4UL , 0 x5UL , 0x2UL , 0 x3UL , 0 x6UL , 0 x7UL } ; __m512i idx_4=(__m512i ) { 0 x8UL , 0 x0UL , 0 x1UL , 0 x2UL , 0x3UL , 0 x4UL , 0 x5UL , 0 x8UL } ; __m512i idx_5=(__m512i ) { 0 x8UL , 0 x8UL , 0 x8UL , 0 x6UL , 0x7UL , 0 x8UL , 0 x8UL , 0 x8UL } ; __m512i idx_6=(__m512i ) { 0 x0UL , 0 x0UL , 0 x4UL , 0 x5UL , 0xcUL , 0 xdUL , 0 x0UL , 0 x0UL } ; __m512i idx_7=(__m512i ) { 0 x0UL , 0 x0UL , 0 x6UL , 0 x7UL , 0xeUL , 0 xfUL , 0 x0UL , 0 x0UL } ;

  m i d d l e = _mm512_permutexvar_epi64 ( idx_2 , R1_512 ) ; m i d d l e ^=_mm512_permutexvar_epi64 ( idx_3 , R2_512 ) ;

  R3_512 ) ; m i d d l e ^= _mm512_permutex2var_epi64 ( R0_512 , idx_7 , R3_512 ) ;

/

  / f i n a l r e c o n s t r u c t i o n ( Karatsuba ) const __m512i perm_cm = ( __m512i ) { 0 x4UL , 0 x5UL , 0x6UL , 0 x7UL , 0 x0UL , 0 x1UL , 0 x2UL , 0 x3UL } ; cm = _mm512_permutexvar_epi64 ( perm_cm , cm ) ; C[ 0 ] = _mm512_mask_xor_epi64 ( c l , 0 x f 0 , c l , cm ) ; C[ 1 ] = _mm512_mask_xor_epi64 ( ch , 0 x 0 f , ch , cm ) ; }

  1)/(x 2 + x) • C(0) + C(1) + C(x)/x +C(x + 1)/(x + 1) + (x 2 + x) • C(∞) c 2 = C(1)/(x 2 + x) + C(x)/(x + 1) + C(x + 1)/x +(x 2 + x + 1) • C(∞) c 3 = C(0)/(x 2 + x) + C(1)/(x 2 + x) + C(x)/(x 2 + x) +C(x + 1)/(x 2 + x) c 4 = C(∞)Appendix E Experimentation procedureMeasurements were performed on a Dell Inspiron laptop with an Intel Tiger Lake processor. vendor_id : GenuineIntel cpu family : 6 model : 140 model name : 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz

  and our work.

	mult.	#pclmul	#xor
	Schoolbook	t 2	4 × t 2
	KaratRec	t log 2 (3)

Table 2 :

 2 Instruction count for the mul128x4 bit multiplication versions

		Instruction count of mul128x4
		128 bit size operands
	Instructions	Drucker et al. [11] schoolbook (this work)
	_mm512_clmulepi64_epi128	3	4
	XOR		4	1
	_mm512_mask_xor_epi64	2	2
	_mm512_permutex_epi64	2	-
	_mm512_permutexvar_epi64	1	1
	Total		12	8
			Instruction count
		512 bit size operands
	Instructions	karat_mult_1_512_SB karat_mult_1_512 SB512 mul512 [11]
	_mm512_clmulepi64_epi128		12	9	16	9
	XOR		11	20	5	22
	_mm512_mask_xor_epi64		11	11	14	11
	_mm512_permutex_epi64		-	-	-	6
	_mm512_permutexvar_epi64		13	20	18	8
	_mm512_permutex2var_epi64		3	3	5	5
	_mm512_alignr_epi64		-	-	-	6
	Total		50	63	58	67

Table 3 :

 3 Instruction count for the 512 bit multiplication versions2.2.8 PerformancesWe now check the performances of all the previous approaches. The test procedure is presented Appendix E page 19. Since the SB512 version is the slowest one, we chose to provide its performances Table14page 20, in Appendix F.

Table 4

 4

	shows that the instruction

Table 4 :

 4 Performance comparison for Algorithm 1 numbers are 24% to 33% lower. This may be explained by the intensive use of the _mm512_clmulepi64_epi128 we have in our case. Due to this instruction high latency and throughput, the lesser instruction count does not yield such a large decrease for the clock cycle numbers. These results may also vary versus the processor version, according to the corresponding _mm512_clmulepi64_epi128 instruction latency and throughput.

Table 5 :

 5 Performance comparison for Algorithm 2 and 3

	, our implementations achieve

Table 6 :

 6 Performance comparison for recursive Algorithms 2 and 3

		size	gf2x	# clock cycles AVX2 AVX512		size	gf2x	# clock cycles AVX2 AVX512
	Karat3(Karat3) Karat5(Karat5)	4608 9216 18432 36864 73728 12800 25600 51200 102400	4129 11261 34205 97586 270900 19495 56258 161740 438283	2707 8102 24498 76912 232362 17988 53132 159148 479256	1920 5425 16340 51097 154065 11471 36417 107020 321361	Karat3(Karat5) Karat5(Karat3)	7680 15360 30720 61440 122880 7680 15360 30720 61440 122880	8317 25110 73222 209872 618612 8242 26364 72897 216280 619802	7075 20670 63386 192557 597357 7031 20586 63653 190317 580575	4496 13845 42454 127336 388414 4943 13775 42713 129419 387763

Table 7 :

 7 762 • 64 -1 = 24191. Complexity of the Toom-Cook's multiplication variants Example 3. We now will use the previous multiplication in order to build a 1-recursive Toom-Cook multiplication, setting w = 256. One now chooses n such that 64n + 2 • 256 = 24192. This leads to n = 370 and t = 1110. This multiplication proceeds polynomials of degree at most 71039:

	mult.	Complexity	size range (gf2x [18])
			operand degrees	# 64 bit words
	KaratRec	O(t log 2 (3) )	< 1343	< 21
	Toom-Cook 3-5	O(t log 3 (5) )	1343 < degrees < 22143	21< w <346
	Toom-Cook 4-7	O(t log 4 (7) )	> 22143	> 346
	Toom-Cook 5-9	O(t log 5 (9) )	above	above

Table 8

 8

	):
	-with KaratRec elementary multiplications whose size
	is among 512, 1024, 2048, 4096, 8192, 16384 and
	32768 bits;

Table 8 :

 8 Performance comparison between AVX2 and AVX512 Toom3 multiplications, with various elementary Karat-

					size	gf2x	# clock cycles AVX2 AVX512
				24192	47200	33501	22238
	Toom3Mult(KaratRec)	48768 144960 104257	69328
				97920 427101 323719 204596
				18048	34019	22032	14872
	Toom3Mult(Karat3)	36480	94378	67451	46940
				73344 272891 206601 137791
				14976	25146	18570	12664
	Toom3Mult(Karat5)	30336 61056 211790 170316 113783 77548 57339 38668
			122496 622830 513595 339387
	suba multiplications
								# clock cycles
	multiplication size				AVX2	AVX512
	Elt. Karat. size Karat3 Toom3 Karat3	Toom3	Karat3	Toom3
	2048	6144	5760		4150		3987	2655	-
	3072	9216	8832		8405		7726	5425	-
	4096	12288	11998	12508	11480	8118	-
	6144	18432	18048	25264	22032	16340	14872
	8192	24576	24192	37024	33501	24488	22238
	12288	36864	36480	79727	67451	51097	46940
	16384	49152	48768	117707 104257	76147	69328
	24576	73728	73344	241548 206601	154065 137791
	32768	98304	97920	357577 323719	231686 204596

Table 9 :

 9 Performance comparison between AVX2 and AVX512 Toom3Mult and Karat3 multiplications

  6: R 1 ← KaratRec(A 1 , B 1 , t/5) 7: R 2 ← KaratRec(A 2 , B 2 , t/5) 8: R 3 ← KaratRec(A 3 , B 3 , t/5) 9: R 4 ← KaratRec(A 4 , B 4 , t/5) 10: R 01 ← KaratRec(A 0 + A 1 , B 0 + B 1 , t/5) 11: R 02 ← KaratRec(A 0 + A 2 , B 0 + B 2 , t/5) 12: R 03 ← KaratRec(A 0 + A 3 , B 0 + B 3 , t/5) 13: R 04 ← KaratRec(A 0 + A 4 , B 0 + B 4 , t/5) 14: R 12 ← KaratRec(A 1 + A 2 , B 1 + B 2 , t/5) 15: R 13 ← KaratRec(A 1 + A 3 , B 1 + B 3 , t/5) 16: R 14 ← KaratRec(A 1 + A 4 , B 1 + B 4 , t/5) 17: R 23 ← KaratRec(A 2 + A 3 , B 2 + B 3 , t/5) 18: R 24 ← KaratRec(A 2 + A 4 , B 2 + B 4 , t/5) 19: R 34 ← KaratRec(A 3 + A 4 , B 3 + B 4 , t/5) 20: // Reconstruction 21: R ← R 0

Table 13 :

 13 Instruction count for the 256 bit multiplication versions

	Instruction count 256 bit size operands		VPCLMULQDQ			AVX2
		Instructions	SB version Karat. version Drucker et al. [10] Karat. Rec.
	_mm512_clmulepi64_epi128	4	3	8
	_mm_clmulepi64_si128						9
		XOR		5	7	15		25
		AND		0	1
	_mm512_broadcast_i64x4	2	2
	_mm512_permutexvar_epi64	3	6	6
	_mm512_permutex2var_epi64	5	5
	_mm512_alignr_epi64				2
	_mm_loadu_si128						4
	_mm_shuffle_epi32						6
	_mm_setzero_si128						6
	_mm_unpacklo_epi64						3
	_mm_unpackhi_epi64						3
		Total		19	24	31		56
	additional vmovdqa64's		8	9	6		8
		KaratRec		Drucker et al.	AVX2	This work
				[10]		[2]	vpclmulqdq-512
	size		gf2x	4x4 -256 8x8 -512		SB-256 Karat.-256
	1024	# clock cycles # instructions	339 1224	239 326	169 221	183 612	167 276	186 339
	2048	# clock cycles # instructions	998 3892	744 1137	532 736	610 1867	541 908	627 1102
	4096	# clock cycles # instructions	2949 11079	2262 3632	1621 2358	1929 5684	1685 2769	1937 3346
	8192	# clock cycles # instructions	8742 33182	6960 11535	4926 7547	6038 17991	5204 8718	6205 10552
	16384	# clock cycles # instructions 100163 26128	21154 35483	14940 23313	18327 54840	15675 26482	18043 32068
	32768	# clock cycles # instructions 295755 78889	65758 108329	45244 70941	59613 166410	47592 80736	54881 97278
	65536	# clock cycles # instructions 853977 226640	203855 328772	140347 214955	187305 503014	147572 244222	169844 293873
	131072	# clock cycles # instructions 2516857 667900	621942 992919	425430 648533	572984 1515625	446811 736345	511939 885164

Table 14 :

 14 Performance comparison for Algorithm 1

We do not present in detail our variant karat_mult_1_512 based on the mul128x4, however, we refer the reader to their paper[START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF] for its presentation.

This work

NIST release [START_REF] Aguilar-Melchior | Hamming Quasi-Cyclic (HQC)[END_REF] improved AVX2 AVX512 We reproduce here the Karatsuba algorithms:

-Algorithm 1 reproduces the recursive multiplication with a two halves split, from [START_REF] Nègre | Impact of Optimized Field Operations AB, AC and AB + CD in Scalar Multiplication over Binary Elliptic Curve[END_REF]; -Algorithm 2 shows the three parts split corresponding approach; -Algorithm 3 shows the five parts split corresponding approach;

Algorithm 1 KaratRec(A,B,t), from [START_REF] Nègre | Impact of Optimized Field Operations AB, AC and AB + CD in Scalar Multiplication over Binary Elliptic Curve[END_REF] Require: A and B on t = 2 r computer words.

return ( M ult64(A, B) ) 3: else 4:

// Split in two halves of word size t/2.

5:

Algorithm 2 Karat3(A,B,t), from [START_REF] Weimerskirch | Generalizations of the karatsuba algorithm for efficient implementations[END_REF] Require: A and B on t = 3 × 2 r computer words. Ensure: R = A × B 1: // Split in three thirds of word size t/3. We present here our schoolbook AVX512 implementation of the 256 bit operand size multiplication, with comments and explanations.

Appendix C Source code for 512-bit operand size multiplications

We detail now the source code of the karat_mult_1_512_SB procedure.

First, the preamble declares the constant indexes for the _mm512_permutex2var_epi64 and _mm512_permutexvar_epi64 instructions. Next, we compute the registers al, ah, bl, bh, sa and sb so that they contain the split parts for the 256 bit operands, and the corresponding sums for the Karatsuba 256 bit middle multiplication. We compute now the three 256 bit multiplications in order to prepare the 512 bit registers cl, ch and cm containing their results.

// f i r s t m u l t i p l i c a t i o n 256 : AlBl These lines computes four 128 bit operand size multiplications in parallel, using a schoolbook approach. This procedure acts like the mul128x4 procedure of Drucker et al. [START_REF] Drucker | Fast polynomial inversion for post quantum QC-MDPC cryptography[END_REF], which is based on the Karatsuba algorithm 2 .

The 512 bit registers l and h now contains the four elementary 256 bit results. We now compute the same the two remaining 256 bit operand size multiplications: // second m u l t i p l i c a t i o n 256 : AhBh . . . ch ^= m i d d l e ;


F De Vleeschouwer

Á G Le Roux

A Hursthouse

D Weiss

over the web where possible

Selected papers from the 29th SEGH Conference on Environmental Geochemistry and Health

Franc ¸ois De Vleeschouwer This special issue represents a selected collection of studies presented during the conference. It covers a wide range of disciplines and approaches applied to the environmental and health impacts of contaminants.

On behalf of the Society for Environmental Geochemistry and Health, we would like to thank all the authors and reviewers who contributed to this special volume issue.

F

  . De Vleeschouwer (&) Á G. Le Roux INP, UPS, EcoLab (Laboratoire Ecologie Fonctionnelle et Environnement), ENSAT, Universite ´de Toulouse, Avenue de l'Agrobiopole, 31326 Castanet Tolosan, France e-mail: francois.devleeschouwer@ensat.fr

  • Andrew Hursthouse • Gae ¨l Le Roux • Dominik Weiss

	The 29th Conference on Environmental Geochemistry	of Speciation, Environmental Behaviour and Human
	and Health took place in Toulouse (France) from 8 to	Health Impacts).
	12 July 2013. It gathered 160 delegates from 35
	countries around regular topics (Biogeochemistry and
	Health, Biogeochemistry and Ecotoxicology, Spatial
	and Temporal Records of Pollution including Catch-
	ment Studies), as well as special Sessions (Frontiers in
	Mercury Biogeochemistry, Bioaccessibility of Pollu-
	tants in Soils and Vegetables, Arsenic: Current Issues