High performance SIMD modular arithmetic
for polynomial evaluation

We are given a multivariate polynomial $f\in{\mathbb{Z}}[x_{1},\dots,x_{n}]$ and we want to evaluate the variables $x_{3},\dots,x_{n}$ in $f$ at integers modulo a prime $p$. Let ${\mathbb{F}}_{p}$ denotes the finite field of integers modulo a prime $p$. The prime $p$ is chosen so that all integer arithmetic in ${\mathbb{F}}_{p}$ can be done with the hardware arithmetic logic units of the underlying machine. For example, with a machine which supports a 64 bit by 64 bit integer multiply, the application would use $p<2^{64}$. We will pick non-zero integers $\beta_{3},\beta_{4},\dots,\beta_{n}$ uniformly at random from ${\mathbb{F}}_{p}$ and compute a sequence of $T$ evaluations

The values $b_{t}$ are polynomials in ${\mathbb{F}}_{p}[x_{1},x_{2}]$. We call them bivariate images of $f$. For convenience we will use the notation

so that we may write $b_{t}(x_{1},x_{2})=f(x_{1},x_{2},\beta^{t})$.

Before presenting two application examples of such computation, we first emphasize that we evaluate here at powers of $\beta$, not at $T$ different random points in ${\mathbb{F}}_{p}^{n-2}$, since this enables one to save computations. Indeed, when Zippel[15] introduced the first sparse polynomial interpolation algorithm and used his algorithm to compute a gcd of two polynomials $f$ and $h$ in ${\mathbb{Z}}[x_{1},\dots,x_{n}]$, he used random evaluation points. To interpolate a polynomial with $s$ terms his method solves one or more $s\times s$ linear systems modulo $p$. Using Gaussian elimination this does $O(s^{3})$ arithmetic operations in ${\mathbb{F}}_{p}$ and requires space for $O(s^{2})$ elements of ${\mathbb{F}}_{p}$. Then Zippel[16] showed that if powers $\beta^{t}$ are used for the evaluation points, the linear systems are transposed Vandermonde systems which can be solved using only $O(s^{2})$ operations in ${\mathbb{F}}_{p}$ and $O(s)$ space. In Sect. 2.3, we will see that using powers of $\beta$ also reduces the polynomial evaluation cost.

Second, we also emphasize that we evaluate at $\beta^{t}$ for $1\leq t\leq T$, and not for $0\leq t<T$, since the evaluation point $\beta^{0}=(1,1,\dots,1)$ may not be usable. For example, consider the following gcd problem in ${\mathbb{Z}}[x_{1},x_{2},x_{3},x_{4}]$. Let

Since $\gcd(c,d)=1$ we have

But suppose we use $\beta^{0}=(1,1)$. Then since $d(x_{1},x_{2},1,1)=0$ we have

We cannot interpolate $g$ using this image. We say $(1,1)$ is an unlucky evaluation point. Such unlucky evaluation points are avoided with high probability by picking $\beta$ at random from ${\mathbb{F}}_{p}^{n-2}$ and evaluating at $\beta^{t}$ for $t$ starting at $1$.

Such bivariate images of $f$ are needed in modern algorithms of Computer Algebra for factoring polynomials with integer coefficients and for computing gcd of polynomials with integer coefficients. Two examples are presented below.

Let $p$ be a prime and let $f\in{\mathbb{F}}_{p}[x_{1},x_{2},\dots,x_{n}]$. We may write $f$ as

where $a_{i}\in{\mathbb{F}}_{p}$ are non-zero, $d_{i}$ and $e_{i}$ are non-negative integers and $M_{i}$ is a monomial in $x_{3},\dots,x_{n}$. For $\beta=(\beta_{3},\beta_{4},\dots,\beta_{n})\in{\mathbb{F}}_{p}^{n-2}$, we want to compute $T$ partial evaluations

If we let $m_{i}=M_{i}(\beta_{3},\dots,\beta_{n})\in{\mathbb{F}}_{p}$ and $M_{i}(x_{3},\dots,x_{n})=\prod_{k=3}^{n}x_{k}^{d_{ik}}$ then we observe that

Thus

Now we can present the “matrix method”[5] which relies on the powers of $\beta$ to efficiently compute the $T$ bivariate images. First we compute the monomial evaluation $m_{i}$ by evaluating $M_{i}(\beta_{3},\dots,\beta_{n})$. To do this, let $d_{k}=\deg(f,x_{k})$ and let $d=\max_{k=3}^{n}d_{k}$. We pre-compute tables of powers

This takes at most $(n-2)(d-1)$ multiplications. Then, for $1\leq i\leq s$ we compute $m_{i}=M_{i}(\beta)$, using $n-3$ multiplications for each $m_{i}$ thanks to the tables of powers, and thus using $(n-3)s$ multiplications in total. Therefore we can compute the $m_{i}$ with $O(nd+ns)$ multiplications. Computing

needs 1 multiplication for each $a_{i}m_{i}\in{\mathbb{F}}_{p}$, hence $s$ multiplications in total. We can compute the next evaluation

using another $s$ multiplications if we save $a_{i}m_{i}\in{\mathbb{F}}_{p}$ for $1\leq i\leq s$ and multiply them by $m_{i}$. This leads to an algorithm that computes the $T$ evaluations using $O(nd+ns)$ multiplications to compute the $m_{i}$, plus a further $sT$ multiplications to compute the $a_{i}m_{i}^{t}$ for $1\leq t\leq T,~{}1\leq i\leq s$, hence $O(nd+ns+sT)$ multiplications in total. With random points instead of powers of $\beta$, the $T$ evaluations would have required a larger operation count of $O(nT(d+s))$ multiplications in total[5].

One way to see the matrix method is to think of evaluating $f$ at $\beta^{t}$ as the following $t\times s$ matrix-vector multiplication.

In practice, the complete matrix is not explicitly built and we take advantage of the connection between successive rows of the matrix. Let $a=[a_{1},a_{2},\dots,a_{s}]$, $m=[m_{1},m_{2},\dots,m_{s}]$ and $X=[x_{1}^{d_{1}}x_{2}^{e_{1}},x_{1}^{d_{2}}x_{2}^{e_{2}},\dots,x_{1}^{d_{s}}x_{2}^{e_{s}}]$. Let $u\circ v$ denote the Hadamard product of two vectors $u,v\in{\mathbb{F}}_{p}^{s}$, that is $u\circ v=[u_{1}v_{1},u_{2}v_{2},\dots,u_{s}v_{s}]\in{\mathbb{F}}_{p}^{s}.$ Then, viewing $b_{1}(x_{1},x_{2})$ and $b_{2}(x_{1},x_{2})$ as vectors of terms, we have

Finally, for a given $t$ we will have to compute the sum $c_{i,t}$ of $a_{i}m_{i}^{t}$ for all $a_{i}m_{i}^{t}x_{1}^{d_{i}}x_{2}^{e_{i}}$ sharing the same $d_{i}$ and $e_{i}$ values. This sum $c_{i,t}$ is indeed the coefficient of $x_{1}^{d_{i}}x_{2}^{e_{i}}$ in $b_{t}(x_{1},x_{2})$. These “coefficient reductions” are required since in Eq. (1), multiple $M_{i}(x_{3},\dots,x_{n})$ can potentially share the same $d_{i}$ and $e_{i}$ values. If the monomials $x_{1}^{d_{i}}x_{2}^{e_{i}}M_{i}(x_{3},\dots,x_{n})$ in the input polynomial $f$ are sorted in lexicographical order with $x_{1}>x_{2}>x_{3}>\dots>x_{n}$ then the monomials $x_{1}^{d_{i}}x_{2}^{e_{i}}$ will be sorted in $X$ which makes adding up $c_{i,t}$ coefficients of like monomials in $b_{t}(x_{1},x_{2})$ straightforward (with $O(s)$ additions for each evaluation). Doing so, we compute $f(x_{1},x_{2},\beta^{t})$ for $1\leq t\leq T$ with $O(nd+ns+sT)$ multiplications and $O(sT)$ additions.

The resulting algorithm for the compute kernel of the matrix method is detailed in Algorithm 1, where we save the successive $a_{i}m_{i}^{t}$ values in the $a$ vector, and where ${\oplus}_{p}$ and ${\otimes}_{p}$ denote the arithmetic operators modulo $p$ ($c=a{\otimes}_{p}b$ denoting $c\equiv a\times b\pmod{p}$, and $d=a{\oplus}_{p}b$ denoting $d\equiv a+b\pmod{p}$ with $(a,b,c,d)\in{\mathbb{F}}_{p}^{4}$).

In this article, we will use by default the following parameters when measuring the time or performance of our partial modular polynomial evaluation with 64-bit integers: $s=5\times 10^{5}$ terms; $n=6$ variables, hence $4$ evaluated variables; a maximum degree of $d=10$ in each variable; and the number of evaluations $T$ chosen here as $10000$ to have a measurable computation time, but $T$ can be much lower in actual use. These parameters have been chosen to be realistic and to lead to stable and reproducible performance results.

Given a fixed¹¹1The value of $p$ is fixed in this article, but is not a constant (from the programming point of view) in our implementation. Our implementation will indeed be used for multiple $p$ values, which are unknown at compile time. The compiler cannot therefore optimize the code for a specific $p$ value. integer $p>1$, we focus on the efficient computation of $c=a{\otimes}_{p}b$ and $d=a{\oplus}_{p}b$ . We target the algorithm that will offer the best performance: such an algorithm must thus be efficient in scalar mode (i.e. non-SIMD), while being also suitable for vectorization. While ${\oplus}_{p}$ can be implemented with a compare instruction, ${\otimes}_{p}$ requires integer divisions which are expensive operations on current processors[23], and for which no SIMD integer division instruction is available in SSE, AVX or AVX-512. Hence, various alternate algorithms have been designed in order to efficiently compute ${\otimes}_{p}$. We briefly recall the most important ones.

In order to compute $c=a{\otimes}_{p}b$, one can first rely on floating-point arithmetic to compute $q=\lfloor\frac{a\times b}{p}\rfloor$ and then deduce $c=a\times b-q\times p$ (see for example Alverson[24], Baker[25]). This requires conversions to/from floating-point numbers, and the number of bits of the floating-point number mantissa has to be twice as large as the number of bits of $p$ (to hold the product). In order to avoid the conversions between floating point numbers and integers, the floating-point reciprocal $p^{-1}$ can be rescaled and truncated into an integer. The quotient $q$ is hence approximated, and some adjustments enable one to obtain the remainder $c$. This integer-based method is known as the Barrett’s product[26, 27]. Another integer-based approach relies on Montgomery’s reduction[28]: a comparison between the two methods has been done for example by van der Hoeven et al.[14]. An improved version of Barrett’s product with integer only operations has been proposed by Möller and Granlund[29]: Monagan and coworkers[5, 17] use an implementation (written by Roman Pearce) of this latter[29] method (with $p^{-1}$ precomputed) in their original code for polynomial evaluation with 64-bit integers. This offers a 11x performance gain[30] for ${\otimes}_{p}$ with respect to one integer division. However, this implementation relies on 128-bit intermediate results: for SIMD processing on 64-bit elements, this implies that only half of the SIMD lanes will be used, hence leading to twice lower SIMD speedups. One can replace these 128-bit variables with two 64-bit variables (hence using only one 64-bit lane per operation), but this requires extra arithmetic and bit shifting. Moreover, to our knowledge there is no SSE/AVX2/AVX-512 intrinsic which performs multiplications on 64-bit integers and provides either the 128-bit results (similarly to the _mm{,256,512}_mul_epu32 intrinsics on 32-bit integers) or their upper and lower 64-bit parts. This greatly complicates the SIMD programming of the Möller and Granlund[29] algorithm for 64-bit integers.

We therefore focus in this article on the use of floating-point (FP) FMA (fused multiply-add) instructions for floating-point based modular arithmetic. Since the FMA instruction performs two operations ($a*b+c$) with one single final rounding, it can indeed be used to design a fast error-free transformation of the product of two floating-point numbers[31]. Such error-free transformation computes the accurate floating-point result of the product. As described and proved by van der Hoeven et al.[14], this makes it possible to design a modular multiplication with double-precision floating-point numbers, provided that $p$ has at most 50 bits: see Algorithm 2. Intuitively, an error-free transformation (Lines 1 and 2 in Algorithm 2) enables one to compute in twice working precision[31], and hence to precisely handle the multiplication result before reduction modulo $p$. More precisely, $\ell$ stores the rounding error of the product $x*y$ (i.e. $h+\ell$ exactly equals $x*y$). The approximate real quotient $(x*y)/p$ is then computed in $b$ using the pre-computed $u=1/(double)p$, and rounded to an (approximate) integer quotient $c$. A first approximate remainder $d$ is computed using $c*p$, and added to $\ell$ in $g$ in order to take into account the initial rounding error. $g$ is finally corrected so that we exactly have: $g\equiv x\times y\pmod{p}$. We emphasize that all this is achieved with 64-bit floating-point numbers only: no larger variables are required and we can thus benefit from the full SIMD speedup (up to 8x on AVX-512). The corresponding FP-based modular addition algorithm is presented in Algorithm 3.

We also emphasize that the limit on the size of $p$ (at most 50 bits) is not problematic regarding our targeted applications (presented in Sect. 2.2). Indeed, let $f(x_{1},x_{2},\dots,x_{n})=\sum_{i=1}^{s}a_{i}M_{i}(x_{1},x_{2},\dots,x_{n})$ be a polynomial we wish to interpolate. Ben-Or/Tiwari [18] and Zippel [16] both pick $\beta\in{\mathbb{F}}_{p}^{n}$ at random and interpolate $f$ from $f(\beta),f(\beta^{2}),f(\beta^{3}),\dots.$ Both methods require the monomial evaluations to be distinct, that is, $M_{i}(\beta)\neq M_{j}(\beta)$ for $1\leq i<j\leq s$. For this to hold with reasonable probability we require $p>100s^{2}$. For a large value of $s$, say $10^{4}<s<10^{6}$, the requirement $p>100s^{2}$ means 32-bit primes are too small but 50-bit primes are sufficient.

Finally, we stress that relying on FMAs is relevant regarding HPC architectures. Current high-end HPC-oriented Intel CPUs with AVX2 or AVX-512 indeed offer two FMA SIMD units. HPC-oriented GPUs from NVIDIA or AMD (not studied in this article), whose performance strongly depends on SIMD computing, also fully support FMA instructions.

We plan to integrate the SIMD implementations of the FMA-based ${\oplus}_{p}$ and ${\otimes}_{p}$ modular operations in the polynomial evaluation algorithm (see Algorithm 1) on AVX2 or AVX-512 CPUs. For this purpose, there are multiple programming paradigms regarding SIMD computing.

A first possibility is to rely on intrinsics programming. Such low-level programming enables the programmer to reach high performance, but at a non-negligible development cost. This will be our primary programming paradigm, and we will detail the corresponding implementations in Sect. 3.3.

A second possibility is to rely on the compiler to benefit from a higher programming level. Compilers can detect parallel and vectorizable loops and automatically vectorize these loops. But, as further detailed in Sect. 3.5, such automatic vectorization will fail in our polynomial evaluation. Therefore, we consider here a third programming paradigm: compiler directives for SIMD programming. In C/C++ programming, these are pragmas which enable the programmer to indicate (and ensure) that a given loop is parallel: no dependency analysis is then required by the compiler. Such compiler directives are available in the Intel C/C++ Compiler ICC (#pragma simd), and have been standardized in the last versions of OpenMP²²2https://www.openmp.org/ (starting from OpenMP 4.0). We will rely here on OpenMP due to its sustainability, its wide usage in HPC, and its availability in both ICC and GCC (the GNU Compiler Collection). Such high level programming with OpenMP will enable us to avoid writing intrinsics, to have one scalar C code for both AVX2 and AVX512, and to avoid array padding or loop splitting when the iteration number is not a multiple of the SIMD vector size. OpenMP directives will also enable us to overcome the limits of the automatic vectorization for our polynomial evaluation. However, the SIMD code generated by the compiler may differ from the intrinsic code and hence lead to lower performance.

In the rest of Sect. 3, we will thus investigate these two SIMD programming paradigms: SIMD intrinsics and OpenMP SIMD directives. Their performance results will also be detailed and compared.

We start with microbenchmarks, presented in Fig. 1, of the two modular arithmetic operations ${\otimes}_{p}$ and ${\oplus}_{p}$ . Like all following performance tests, these microbenchmarks have been performed on the two compute servers AVX2 server and AVX-512 server presented in Table 1. We first notice that, regarding the original implementation of the modular multiplication (integer based, by R. Pearce), our microbenchmark results are consistent with the 6.016 cycles obtained with GCC by Monagan and coworkers[30] on older CPUs.

Regarding ${\otimes}_{p}$ (see Figs. 1(a) and 1(b)), and with respect to the original integer based implementation, the scalar floating-point (FP) based implementation offers lower performance when including back and forth conversions between integers and floating-points numbers, and similar performance when not considering these conversions. The SIMD FP based implementation offer (with intrinsics, and without considering the conversions) a 3.2x (resp. 3.7x) speedup with GCC (resp. ICC) over the scalar FP-based implementation on AVX2 server. On AVX-512 server, the speedup is 5.9x (resp. 7.2x) speedup with GCC (resp. ICC). This shows that the performance gain of our new AVX-512 ${\otimes}_{p}$ implementations is indeed twice greater than the AVX2 one.

Regarding ${\oplus}_{p}$ (see Figs. 1(c) and 1(d)), the scalar FP-based implementation leads to much greater cycle numbers than the integer-based one: this is due to the branching of the compare instruction required in the FP implementation (see Algorithm 3). The comparison in the integer-based ${\oplus}_{p}$ implementation can indeed be replaced by shifting[34, 14], which avoids the branching performance impact on the pipeline filling. Thanks to the use of the AVX2 blendv_pd intrinsic and of AVX-512 masks, there is also no branching in the SIMD FP-based implementations (with intrinsics) which implies a strong performance gain (around one order of magnitude), in addition to the SIMD speedup. With respect to the scalar integer-based one, the SIMD speedups of the FP-based implementations with intrinsics (without considering the conversions) are 4.0x (resp. 3.6x) with GCC (resp. ICC) on AVX2 server, and 8.7x (resp. 8.5x) with GCC (resp. ICC) on AVX-512 server. This shows that our AVX-512 ${\oplus}_{p}$ implementation is twice faster than the AVX2 one of Van der Hoeven et al.[14].

When considering the conversions, the overhead of these conversions can annihilate the SIMD performance gain on AVX2 server. This is due to the lack of AVX2 conversion instruction between 64-bit integers and 64-bit floating-point numbers: the conversions are thus performed in scalar mode which has a strong performance impact. In comparison, such a SIMD instruction is available in AVX-512³³3More precisely, the 64-bit conversions belong to the AVX-512DQ instruction set which is available on our Intel Xeon Gold 6152 CPUs, but not on the prior Intel Knights Landing (Xeon Phi) processors. , where conversions can be performed in SIMD mode. Their overheard is therefore much lower on AVX-512 server.

Finally, we also consider using OpenMP to vectorize the code. The first issue lies in having the compiler generate SIMD FMA instructions from the fma() function call in the C+OpenMP code for ${\otimes}_{p}$. This is effective with GCC thanks to the -fno-trapping-math option which allows us to assume that floating-point operations cannot generate traps, such as division by zero, overflow, underflow, inexact result and invalid operation. Unfortunately, using all possible floating-point model variations (-fp-model options) did not enable us to generate SIMD FMA instructions with ICC. The ICC OpenMP code hence relies on scalar FMA instructions, which explains its important performance overhead over the GCC OpenMP code for ${\otimes}_{p}$. Secondly, on AVX-512 server we had to force the AVX-512 vectorization using -qopt-zmm-usage=high with ICC and -mprefer-vector-width=512 with GCC, otherwise only AVX2 instructions are generated.

As far as ${\otimes}_{p}$ is concerned, the OpenMP code (with GCC) has the same performance than the SIMD code written in intrinsics on AVX-512 server, but is slower by 18% on AVX2 server. This is because of one additional compare instruction added by the compiler before each blendv_pd instruction. This comparison to zero (either $g>0$, or $g-p\geq 0$) is here to prevent any issue with IEEE 754 signed zeros and the blendv_pd instruction. The compiler is indeed unaware of our specific context which enables us not to consider $-0.0$, as shown in Sect. 3.3.

For ${\oplus}_{p}$, the compiler similarly adds one unnecessary compare instruction which results in a 45% performance penalty on AVX2 server for the OpenMP code (with GCC). However no branching instruction is generated in the SIMD code for ${\oplus}_{p}$ with OpenMP, which makes this OpenMP still rather efficient with respect to the scalar integer-based implementation: 2.8x faster on AVX2 server, and 7.4x on AVX-512 server (with GCC).

We can now consider the integration of SIMD modular arithmetic in our partial polynomial evaluation. Due to the cost of the conversions between integers and floating-points numbers (see Sect. 3.4), we choose to perform the first conversion (from integers to floating-point numbers) for each value of the $a$ and $m$ vectors once before the evaluation (i.e. just before Line 1 in Algorithm 1). These conversions are performed in-place to save memory. The reverse conversion (from floating-point numbers to integers) is only performed once for each reduction result (i.e. the $c$ value at Line 10 in Algorithm 1).

Figure 2 presents performance results for our polynomial evaluation using various modular arithmetic implementations. One can see that the scalar floating-point based modular arithmetic makes our polynomial evaluation about 2.5 times slower than the original implementation by Monagan and coworkers[5] (using integer-based modular arithmetic). This is due to the slow FP-based ${\oplus}_{p}$ implementation (because of its branch instruction: see Sect. 3.4).

We now consider SIMD intrinsics to integrate SIMD FP-based modular arithmetic in our originally scalar polynomial evaluation (see Algorithm 1). The resulting SIMD algorithm is written in Algorithm 6. First, a reduction has to be computed within the SIMD vector at the end of the inner loop (Line 14 in Algorithm 6), in order to obtain the final scalar $c$ value. Instead of performing a sequential reduction with the scalar ${\oplus}_{p}$ and its branch instruction, we use SIMD shuffle instructions to write a parallel tree-shaped reduction using only SIMD ${\oplus}_{p}$ operations. Second, we also have to consider memory alignement which can be important for efficient vector loads and stores. However, the $j$ indices used in the inner loop (Line 7 in Algorithm 6) do not lead to aligned memory accesses since the successive $J$ values are not necessarily multiples of the SIMD width. One could choose to make a copy of the vectors $a$ and $m$ with relevant zero padding to ensure aligned memory accesses, but this would require twice the memory. In order to obtain good SIMD speedups without padding, we explicitly decompose this inner loop into three successive loops $L1$, $L2$ and $L3$ (not shown in Algorithm 6). $L1$ and $L3$ have an iteration count lower than the SIMD width and are vectorized thanks to explicit masks. These two loops ensure that the vectorized $L2$ loop perform aligned memory accesses (i.e. its first $j$ index is a multiple of the SIMD width). We noticed that the use of such a vectorized reduction and of such vectorized $L1$ and $L3$ loops offer additionnal performance gains when processing a lower number of terms: e.g. up to 29% for $s=5\times 10^{4}$ terms. With respect to the original implementation using integer-based arithmetic, the resulting polynomial evaluation with SIMD intrinsics offers performance gains ranging from 3.13x (resp. 3.89x) with GCC (resp. ICC) on AVX2 server to 4.18x (resp. 4.74x) with GCC (resp. ICC) on AVX-512 server (see Fig. 2).

Finally, we also consider using OpenMP vectorization for the SIMD FP-based modular arithmetic in our polynomial evaluation. We rely on the new declare reduction directive (available since OpenMP 4.0) to instruct the compiler that the final reduction (Line 14 in Algorithm 6) has to be performed using our modular arithmetic. We emphasize that the vectorization is achieved here thanks to this OpenMP directive (along with the directive which instructs the compiler to vectorize the loop with a reduction). Without such directives given by the programmer, the GCC and ICC compilers both manage to vectorize the microbenchmarks presented in Sect. 3.4 (with the same performance as the OpenMP version), but both fail to vectorize the polynomial evaluation code (because of the required specific reductions). As shown in Fig. 2, the ICC OpenMP vectorization is inefficient due to the FMA issue (see Sect. 3.4), whereas the GCC OpenMP vectorization offers computation times somewhat slower than the intrinsic vectorization: 22% slower on AVX2 server, and 8% on AVX-512 server. It can also be noticed that GCC currently fails to generate align memory accesses (despite the use of the aligned OpenMP clause) and SIMD ${\oplus}_{p}$ reductions, as we do with intrinsics.

In conclusion, while the use of scalar FP-based modular arithmetic lowers the performance of the polynomial evaluation, the SIMD FP-based modular arithmetic clearly improves its performance (up to 4.74x). In the rest of the article, we will rely on the SIMD implementation with intrinsics, and not on the OpenMP one. This is due to the OpenMP performance issue with ICC and to the somewhat lower performance of the SIMD code generated with OpenMP, especially on AVX2 which still equips the vast majority of available CPUs at the time of writing. We however emphasize that the performance results of OpenMP with GCC on AVX-512 are very promising for the future and show the relevance of this approach.

As a last remark, we recall that regarding the microbenchmarks presented in Sect. 3.3 the AVX-512 ${\otimes}_{p}$ and ${\oplus}_{p}$ implementations are twice as fast than the AVX2 ones. Here however for the polynomial evaluation, the AVX-512 performance is only 1.57x (resp. 1.47x) faster than the AVX2 one with GCC (resp. ICC). We believe that this is due to the difference in operational intensities. Indeed, the microbenchmarks performed in Sect. 3.3 have been intentionally designed to be compute-bound in order to measure the number of cycles of the arithmetic operations, and not of the memory accesses. But the operational intensity of our polynomial evaluation is much lower: the Hadamard product and the coefficient reduction correspond to a dot product which is a memory-bound operation in classic floating-point arithmetic. More precisely, the floating-point based modular arithmetic requires 9 flop (floating-point operation) for ${\otimes}_{p}$ (see Algorithm 2) and 2 flop for ${\oplus}_{p}$ (see Algorithm 3), versus 3 memory accesses for each (2 loads and 1 store, without considering $u$ and $p$). This makes our compute kernel (i.e. our polynomial evaluation) less memory-bound than a floating-point dot-product, but the operational intensity of our kernel is not high enough to make it compute-bound: memory accesses are still important in the kernel performance. These memory accesses also tend to lower the performance gain due to the increased compute power of the AVX-512 SIMD units, with respect to the AVX2 units, since there is more stress on memory bandwidth with AVX-512 instructions than with AVX2 ones. We will show how to increase this operational intensity and the polynomial evaluation performance in the next section.

We first rely on the consecutive powers of $\beta$ used for the successive evaluations in the matrix method (see Sect. 2.3). Hence in Algorithm 6, if we consider two consecutive polynomial evaluations $t$ and $t+1$, the values computed in the $a$ vector for the evaluation $t$ are re-used as input for evaluation $t+1$. But for large $s$ values ($s$ denoting the number of terms, see Sect. 2.3), the $a$ elements may have been moved out of the CPU caches. We hence consider computing multiple evaluations at a time, and we denote by $T_{d}$ the number of such (dependent) evaluations. $T_{d}$ is an algorithmic constant, known at compile time. We can then explicitly avoid storing and reloading data from the vector $a$ to/from memory between these $T_{d}$ evaluations. We can also load only once $m$ data from memory for these $T_{d}$ evaluations. Such data reuse increases the operational intensity of our kernel by reducing the number of memory accesses.

However each evaluation depends on the output of the previous one. Even if some operations can be performed concurrently (such as the ${\oplus}_{p}$ operation of the $t$ evaluation and the ${\otimes}_{p}$ of the $t+1$ evaluation), this dependency limits the instruction-level parallelism, and hence prevents us from filling the pipelines of the floating-point units.

Therefore, we rewrite the loop over the $T$ evaluations (Line 1 in Algorithm 1) in order to have fully independent polynomial evaluations. For this purpose, we adapt the algorithm used by Hu and Monagan[5] to introduce thread-level parallelism on multi-core CPUs (see Sect. 2.3.1) in order to introduce here instruction-level parallelism in our compute kernel. Namely, denoting by $T_{i}$ the desired number of independent evaluations (like $T_{d}$, $T_{i}$ is an algorithmic constant, known at compile time), we first precompute $\Gamma=[m_{1}^{T_{i}},m_{2}^{T_{i}},\dots,m_{s}^{T_{i}}]$ using $O(s\log_{2}{T_{i}})$ SIMD multiplications. We also precompute ${\Lambda}_{k}=a\circ[m_{1}^{k+1},m_{2}^{k+1},\dots,m_{s}^{k+1}]$ for $0\leq k<{T_{i}}$ using $O(sT_{i})$ SIMD multiplications. Then, for the computation of the $T$ evaluations we will first perform the coefficient reductions for the first $T_{i}$ evaluations (i.e. on $({\Lambda}_{k})_{0\leq k<{T_{i}}}$), then the Hadamard product ${\Lambda}_{k}\leftarrow{\Lambda}_{k}\circ\Gamma$ (with $0\leq k<{T_{i}}$) for the second chunk of $T_{i}$ evaluations. This will be repeated (coefficient reductions on the previous $T_{i}$ evaluations, then Hadamard product for the next $T_{i}$ evaluations) until all $T$ evaluations have been processed.

This instruction-level parallelism helps fill the instruction pipelines. Moreover the ${\oplus}_{p}$ and ${\otimes}_{p}$ operations are now inverted. Contrary to Algorithm 1 (Lines 7-8) where the second operation depends on the output of the first one, the second operation now only depends on the input of the first one. The two operations can thus be more overlapped, hence easing the pipeline filling.

The main drawback of using $T_{i}$ independent evaluations is the extra memory requirements. For each independent evaluation $k$ we have to store an extra copy $\Lambda_{k}$ of the complete $a$ vector. Moreover, for each independent evaluation $k$ we have to load the $\Lambda_{k}$ vector from memory and store its update in memory. The operational intensity is thus only improved for the $\Gamma$ memory accesses. Therefore, introducing an extra independent evaluation increases less the operational intensity than introducing an extra dependent evaluation.

There is hence a trade-off between pipeline filling and operational intensity regarding the numbers of dependent ($T_{d}$) and independent ($T_{i}$) evaluations. We will thus consider an algorithm where we introduce $T_{d}$ dependent evaluations for each of the $T_{i}$ evaluations, hence computing together $T_{i}\times T_{d}$ evaluations at a time. The loop over $T_{d}$ is chosen as the outer one, and the loop over $T_{i}$ as the inner one: this results in better performance than the opposite loop ordering, which indicates that pipeline filling is here more important than increasing the operational intensity. The optimal values for $T_{d}$ and $T_{i}$ depend on the compiler and on the CPU hardware features, and these will have to be determined in practice using parameter testing and tuning. Since the total number of evaluations $T$ is not necessarily a multiple of $T_{i}\times T_{d}$, the remaining evaluations are processed first by blocks of $T_{i}\times N^{\prime}_{d}$ evaluations (with $1\leq N^{\prime}_{d}<T_{d}$) and then with $N^{\prime\prime}_{d}$ dependent evaluations (with $N^{\prime\prime}_{d}<T_{i}$), where  $T\,mod\,(T_{i}\times T_{d})=T_{i}\times N^{\prime}_{d}+N^{\prime\prime}_{d}$ .