Alexandre Reiffers-Masson
email: alexandre.reiffers-masson@imt-atlantique.fr

Rbccps

Worst-case drift detection of sensor networks: performances and algorithms

The goal of this paper is to derive algorithms that are able to detect unreliable/drifted sensors, in a sensor network. To recover the state of each sensor, we restrict ourselves to a particular decoder, inspired by graph partitioning problems. We provide necessary and sufficient conditions over the measurements such that the decoder perfectly recovers each sensor binary state. The outputs of the decoder can be computed using a dynamic programming approach. One challenging part of this approach is the complexity of the dynamic programming equation. Indeed, the resolution time will increase exponentially with the number of sensors. Therefore, we propose an efficient heuristic method that approximately solves the problem. We study the performance of our algorithm using simulations.

I. INTRODUCTION

In a low-cost sensor network, sensors can gradually become unreliable. This phenomenon is called drift. For instance, it can arise either due to the loss of sensitivity (ex: electrochemical sensors) or aging of the material after long exposure to meteorological conditions. Such phenomena are observed in low-cost pollution sensors. In this paper, the task that we are trying to solve is the detection of drifted sensors when only the measurements of the sensor network are available.

Multiple papers have been trying to solve this problem and we describe the most relevant ones, due to the lack of space. Most drift estimation techniques depend on probabilistic filtering methods. For example in [START_REF] Ms Takruri | Recursive bayesian approaches for auto calibration in drift aware wireless sensor networks[END_REF] and [START_REF] Kumar | Automatic sensor drift detection and correction using spatial kriging and kalman filtering[END_REF] the authors propose a method that interpolates ground truth based on spatial interpolation techniques such as kriging. Drifts were estimated based on the differences in the measured sensor readings and the kriging prediction, which were filtered using April 23, 2020 DRAFT a Kalman filter. However, the accuracy of these two methods is limited by the fact that they accumulate errors over time. We take a different approach, and solve a clustering problem by exploiting the band limited nature of the sensing field, as in [START_REF] Wang | Blind drift calibration of sensor networks using sparse bayesian learning[END_REF], and we treat degraded sensors as outliers in our clusters as in [START_REF] Katz | Outlier-robust convex segmentation[END_REF] and [START_REF] Toby | Clusterpath an algorithm for clustering using convex fusion penalties[END_REF].

In this paper, our contributions are the following: First, we introduce a particular decoder designed to detect drifted sensors. This decoder is the solution of a temporal graph partitioning problem. Secondly, we provide necessary and sufficient conditions over the measurements such that our decoder perfectly recovers the state of each sensor. Thirdly, we propose an efficient heuristic method that approximately solves the temporal graph partitioning problem.

II. SETUP

1) Time series model: Consider a network of I sensors that measures a finite set of spatially and temporally varying quantities. The set of sensors is denoted by I := {1, . . . , I}. Let x n i ∈ R be the output of the sensor i ∈ I, at a particular instant n ∈ {1, . . . , N }. We assume that there are C underlying processes. We denote by s n c one of such process, with c ∈ C := {1, . . . , C}. At each instant n, we consider a scenario where the set of sensors can be partitioned into two sets

I 1 (n) and I 2 (n) such that I 1 (n) ∩ I 2 (n) = ∅ and I 1 (n) ∪ I 2 (n) = I. At each instant n, I 1 (n) is
the set of reliable sensors and I 2 (n) is the set of drifted sensors. We assume that once a sensor drifts, it is unreliable for the remaining time. The latter can be translated into the following condition, I 2 (n) ⊆ I 2 (n + 1). At every instant n, the number of unreliable sensors is smaller that the number of reliable sensors, i.e. | I 1 (n) |>| I 2 (n) |, where | A | denote the cardinality of a set A. We assume that a sensor measures only one underlying process and we consider a linear dynamical model as shown below:

x n i = s n c(i) + n i , if i ∈ I 1 (n), x n i = s n c(i) + n i + δ n i , if i ∈ I 2 (n), (1)
where for each i ∈ I, c(i) ∈ C denote the underlying process measured by sensor i, n i ∈ R denote the noise of sensor i at instant n. For i ∈ I 2 (n), δ n i ∈ R captures the drift value for sensor i at time n.

2) Noise and drift models: At every instant n and for every sensor i, we assume that the noise is such that

( n i ) 2 ≤ c . Moreover, we assume that if i ∈ I 2 (n), then c δ ≤ (δ n i ) 2 ≤ c δ .
We restrict our attention to deterministic approaches for the performance analysis. Indeed, we want to understand the performance of the drift detection task in worst-case scenarios. More precisely, [START_REF] Wang | Blind drift calibration of sensor networks using sparse bayesian learning[END_REF] DRAFT what are the conditions on c , c δ and c δ such that the drifts can be detected and the time series can be filtered/recovered. Such approaches have been proven to be related to some specific noises model (Gaussian noise models for instance, see [START_REF] Başar | H-infinity optimal control and related minimax design problems: a dynamic game approach[END_REF]).

3) Side information: Throughout this paper, we assume that, for every (c, i) ∈ C × I, the underlying process c measured by sensor i is unknown (therefore for every i, c(i) is also unknown).

However, we assume that it is possible to access a similarity matrix E := [[e ij ]] 1≤i,j≤I which captures the similarity between sensors. For instance, if past reliable measurements x -N:0 := [[x n i ]] 1≤i,j≤I,-N ≤n≤0 are available, in the sense that for every n ∈ {-N, . . . , 0}, I 2 (n) = ∅, then it is possible to construct E using the Radial Basis Function (e ij = exp(-

x -N:0 i -x -N:0 j 2 /σ)).
If the similarity matrix is able to recover the entire underlying processes, then e ij = 1 c(i)=c(j) . We will assume a less restrictive assumption, which is

E = E * + M E where; (1) E * = [[e * ij ]] 1≤i,j≤I with e * ij = 1 c(i)=c(j) , (2) M E ∈ M I×I (R), such that max ij {M Eij } ≤ c M E .
We assume that E is a symmetric matrix.

III. DRIFTS DETECTION: APPROACH AND FEASIBILITY

A. Approach

The task studied in this paper is the following: How to recover p n i := 1 i∈I 1 (n) for every sensor i and every instant n knowing x 1:N := [[x n i ]] 1≤i,j≤I,1≤n≤N and E? To solve such challenge, we reformulate our task as the problem of recovering binary sensors states from noisy edge measurements in a temporal graph. The definition of our temporal graph is as follows: the set of sensors I are the nodes of the graph; and an edge at instant n between sensor i and sensor j is given by w n ij = (x n i -x n j ) 2 . Depending on which set (I 1 (n) or I 2 (n)) the node i and the node j are belonging to, the edge w n ij will have different values and therefore will reveal different information about the state of node i and node j:

w n ij =                                ( n i -n j ) 2 , if (i, j) ∈ (I 1 (n)) 2 , ( n i -n j ) 2 + (δ n i ) 2 +2δ n i ( n i -n j ), if i ∈ I 2 (n), j ∈ I 1 (n), ( n i -n j ) 2 + (δ n i -δ n j ) 2 +2(δ n i -δ n j )( n i -n j ), if (i, j) ∈ (I 2 (n)) 2 .
(

) 2
April 23, 2020 DRAFT

Having observed a set of matrices {W n : 1 ≤ n ≤ N } and an adjacency matrix E, the receiver uses a decoder D : (M I×I (R + )) N × M I×I (R + ) → {1, . . . , N } I , where the decoder

D({W n : 1 ≤ n ≤ N }, E) is returning a vector m := [m i ] 1≤i≤I .
For each sensor i, m i + 1 is the predicted time at which sensor i starts to drift. We say that the decoder D is able to detect the drifts in the sensor network if

D({W n : 1 ≤ n ≤ N }, E) = [ N n=1 1 p n i =1 ] 1≤i≤I =: m * .
In this paper, we will restrict to a particular decoder, inspired from graph partitioning problems.

Consider the decoder D λ such that D λ ({W n : 1 ≤ n ≤ N }, E) is the optimal solution of the following integer programming problem:

min m∈{1,...,N } I C(m) := N n=1 I i=1 I j=1 e ij 1 m i ≥n 1 m j ≥n (w n ij -λ). (3)
One important observation concerning the optimization problem ( 3) is that each of its local minimums (denoted by m loc ) will satisfy the following local optimality condition:

m loc i n=1 I j=1 e ij 1 m loc j ≥n (w n ij -λ) ≤ m n=1 I j=1 e ij 1 m loc j ≥n (w n ij -λ), ∀i ∈ I, ∀m = m loc i . (4)
The term e ji 1 m loc j ≥n (w n ji -λ) does not appear in ( 4), because the matrix E and the matrix W n are symmetric matrices.

Our decoder has an elegant interpretation in terms of graph partitioning. For every instant n and for a sensor set J n ⊆ I, we define the weight of the set J n by Ω n (J n ) = (i,j)∈J n e ij (w n ij -λ). The optimization problem defined in (3) is equivalent to:

min {J n ⊆I; n∈{1,...,N }} N n=1 Ω n (J n ), s.t. J n+1 ⊆ J n , ∀n < N.
This optimization problem is about finding the sequence of sets of sensors J n such that the sum of weights e ij (w n ij -λ) inside J n is as small as possible. Therefore we expect the decoder to decide that the sensors with the smallest w n ij will be inside the set J n . The constraint J n+1 ⊆ J n ensures that once a sensor is removed from the set J n , it should be removed forever. The goal of this constraint is to introduce robustness and to ensure that the decoder is considering a sensor as an unreliable sensor at time n if in the future (n > n) the decoder will also take the same decision.

B. Going back to the time series model

In the rest of the paper, we will derive assumptions and algorithms for the general problem of recovering binary sensors states from noisy edge measurements in a temporal graph. But as an April 23, 2020 DRAFT exercise, we show how easy it is to go from assumptions over the matrices W n to assumptions over the time series.

If we are going back to the time series model described in section II, our approach can be understood as follows. Let us assume that c δ -2 √ c c δ > 2c . Then by taking c δ -2 √ c c δ > λ > 2c , we will have w n ij -λ < 0, for every (i, j) ∈ (I 1 (n)) 2 and w n ij -λ > 0, for every (i, j) ∈ I 1 (n) × I 2 (n). In this case the only global optimum of ( 3) is m * . But we believe that this assumption is too strong and in the next section, we derive necessary and sufficient conditions for perfect recovery when it exists (i, j)

∈ (I 1 (n)) 2 such that w n ij > λ and (i, j) ∈ I 1 (n) × I 2 (n) such that w n ij > λ.
For the rest of the paper, all our results will be over the matrices W n , but the reader should keep in mind that it is always possible to go back to the time series model.

C. Main Results

This part of our paper establishes necessary and sufficient conditions for the decoder (3) to be able to perfectly recover the state of each sensor. We define the set

C(i) := {j ∈ I | c(j) = c(i)}, the set A(i, n) := {j ∈ I -{i} | j ∈ I 1 (n) and c(j) = c(i)} and the set B(i, n) := {j ∈ I -{i} | j ∈ I 2 (n) and c(j) = c(i)}.
We assume the following:

Assumption A 1) If i ∈ I 1 (n), then (| A(i, n) |) -1 j∈A(i,n) w n ij ≤ Φ 1 . 2) If i ∈ I 2 (n), then Φ 2 ≥ (| A(i, n) |) -1 j∈A(i,n) w n ij ≥ Φ 3 .
3) For every (i, j)

∈ (I 1 (n)) 2 , w n ij ∈ [c 1 , c 1 ], for every (i, j) ∈ I 1 (n) × I 2 (n), w n ij ∈ [c 2 , c 2 ] and for every (i, j) ∈ (I 2 (n)) 2 , w n ij ∈ [c 3 , c 3 ]. Moreover c 1 ≤ c 2 ≤ c 1 ≤ c 2 , c 1 ≤ c 3 ≤ c 1 ≤ c 3 .
It is important to observe that Assumption A is only used in the analysis of the decoder, it is not needed for the implementation of any algorithm solving (3).

Theorem 1. Necessary conditions:

If assumption A is satisfied, and if

         (1 -α)(Φ 1 -λ) | A(i, n * 1 ) | +αC M E (max{c 2 , c 3 } -λ) | C(i) | +αC M E (c 2 -λ) j =i | C(j) |≤ 0, (1 -α)(Φ 3 -λ) | A(i, n * 1 ) | +αC M E (c 2 -λ) | C(i) | +αC M E (c 2 -λ) j =i | C(j) |≥ 0, (5)
with n * 1 = argmax n A(i, n), then m * i is a local optimum of (4).

April 23, 2020 DRAFT

Sufficient conditions:

If assumption A is satisfied, if (5) is also satisfied and finally if

                   (1 -α) ((Φ 1 -λ) | A(i, n * 2 ) | +(c 2 -λ) | B(i, n * 2 ) |) +αC M E (max{c 2 , c 3 } -λ) | C(i) | +αC M E (max{c 2 , c 3 } -λ) j =i | C(j) |≤ 0, (1 -α) ((Φ 2 -λ) | A(i, n * 3 ) | +(c 3 -λ) | B(i, n * 3 ) |) +αC M E (max{c 2 , c 3 } -λ) | C(i) | +αC M E (max{c 2 , c 3 } -λ) j =i | C(j) |≤ 0, (6)
where

n * 2 = argmax n (Φ 1 -λ) | A(i, n) | +(c 2 -λ) | B(i, n) | and n * 3 = argmax n (Φ 2 -λ) | A(i, n) | +(c 3 -λ) | B(i, n) |, then m *
i is the unique global optimum of (3).

Proof. See appendix.

IV. ALGORITHMS

In this section, we introduce two algorithms. The first one finds the global solution of (3), however, it is not useful if I > 10. The second algorithm finds a local solution of (3), but is much faster.

A. Global solution: Dynamic Programming

The decoder (3) is a nonlinear integer programming problem and it is interesting to note that this problem can be recast as a finite horizon Markov decision problem as stated in the next theorem.

Theorem 2. Let p ∈ {0, 1} I and n ∈ {1, . . . , N }. The solution of the decoder (3) is also the solution of the following dynamic programming equation, ∀n ∈ {1, . . . , N }:

         V (p, n, 1 p < 1) = min u i ∈{0,1}, ∀i (i,j)∈I 2 e ij +V (up, n + 1), ∀p ∈ {0, 1} I , (7)
with the final condition V (p, N + 1) = 0 for all p ∈ {0, 1} I , with e ij = e ij (w ij -λ)u i u j p i p j and starting from the state (1, 1).

Proof. We first introduce the new variable u n i ∈ {0, 1} for all i ∈ I and n ∈ {1, . . . , N }. Then m i = N n=0 p n i with:

p n i = u n i p n i , ∀i ∈ I, ∀n ∈ {1, . . . , N }, p 0 i = 1, ∀i ∈ I.
April 23, 2020 DRAFT Therefore we have a deterministic system where the state is denoted by p n i and the control given by u n i . By simply looking at the cost function of (3), the statement of the theorem follows.

The previous theorem seems to indicate that we could use all the reinforcement learning techniques to solve the drift detection problem. However, it is not feasible in a real-world setting because the number of states in the dynamic programming equation is equal to N × 2 I , and most of the time, we will face situations where I > 10. Therefore, in the next section, we will derive a new algorithm, to obtain an approximate solution of (3).

B. Local solution: Ordinal Potential game

Knowing the complexity of solving the dynamic programming equation ( 7), we provide an algorithm which converges to a local optimum of (3). When the sufficient conditions of theorem 1 are satisfied, the proposed algorithm will converge to the global optimum. Our algorithm is based on a game theory reformulation of (3). The intuition is described below. We assume that every sensor/agent i is minimizing the following local function:

U i (m i , m -i ) = m i n=1 I j =i e ij 1 m * j >0 (w n ij -λ).
A Nash equilibrium of this game is defined as the vector m N E := [m N E i ] 1≤i≤I where for every i, m N E i satisfies:

U i (m N E i , m N E -i ) ≤ U i (m, m N E -i ), ∀m = m N E i ,
where m -i = [m j ] j∈I-{i} . We are now able to prove that the previous optimization problem [START_REF] Katz | Outlier-robust convex segmentation[END_REF] can be seen as an ordinal potential (see [START_REF] Monderer | Potential games[END_REF] for a better understanding of potential games) of our game:

U i (m i , m -i ) > U i (m i , m -i ) ⇔ C(m i , m -i ) > C(m i , m -i ).
Therefore two observations can be made: (1) Each Nash Equilibrium is a local optimum of (3);

(2) One-sided better reply dynamic will converge to a Nash equilibrium. The one-sided better reply dynamic is the following:

One-sided better reply algorithm 1) Set T .

2) For t = 1, . . . , T do:

April 23, 2020 DRAFT a) For i = 1, . . . , I do:

m i (t) = argmin m i ∈{1,...,N } U i (m i , m -i (t)), b) Set m(t + 1) = m(t).
3) Return m(T ).

We do not provide a proof for the number of iterations needed for the convergence of the one-sided better reply dynamic, but we observe that in all our simulations (see section V), the convergence of the algorithm happens in less than 1s.

V. NUMERICAL STUDIES

This section is dedicated to the numerical studies performed to understand the robustness of our approach. We will first describe our two setups. Then we will discuss the different results obtained. Note that, in both scenario we assume that M Eij ∼ N (0, 0.01) and M Eij = M Eji , for all i, j. As a rule of thumb, we choose λ to be equal to the sample quantile, corresponding to the probability 0.7, of the vector [w n ij ] i,j,n . Through different simulations, this approximation has been proved to be efficient. In subsequent work, we will discuss what are the different options for λ. In all the numerical studies, we have only used the one-sided better reply algorithm.

A. Description of the simulation setups

Setup (1): In the first setup, we assume that C := {1, 2} and s n 1 ∼ N (0, 0.1) and s n 2 ∼ N (1.2, 0.1) for every n. Moreover c(i) = 1 for all i ∈ {1, . . . , I/2 } and c(i) = 2 otherwise.

Concerning the drift, sensors are picked randomly as well as the instant at which a sensor starts to drift. Finally if the sensor i drifts, then x n i ∼ N (µ i , 0.1) for every n. In this setup, we choose µ i = (2x i -1)(i + 1), where x i ∼ Bernouilli(0.5). The major purpose of this setup is to prove and to ensure the performance of our algorithm. It does not have any ambition of capturing a real-world scenario. An illustration of this setup, with four drifted sensors, is presented in figure 1.

Setup (2):

In the second setup, we will try to capture a more realistic scenario. Here, we assume that C := {1, 2}, the source signals are ARMA processes for which the parameters are also chosen randomly. As in [START_REF] Wang | Blind drift calibration of sensor networks using sparse bayesian learning[END_REF], we assume that the classical low-cost sensor drift model is similar to the increase in the variance due to the gradual acquisition of drift. We introduce degradation at a randomly selected time n for the sensor i as δ n i = δ n i + n i , where n i ∼ N (0.01, 0.01). Again, the drift sensors are picked randomly. This setup with four sensors, is illustrated in figure 2.

B. Results

We define two performance metrics: The average detection delay and the classical rate of False Positive/Negative. If a sensors are designed to fail in a predetermined number of trials, the average delay is the mean of the absolute difference between the true instant of drift acquisition and the estimated time of degradation. Fig. 3 represents the results of the One-sided better reply algorithm. For each point on the x-axis, the detection time and the rate of False Positive/Negative are averaged over 10 trials. As expected, the decoder is performing better in setup 1 than setup 2.

April 23, 2020 DRAFT

In setup 1, we observe an unexpected phenomenon. The rate of false-negative decreases when the number of unreliable sensors increases. Otherwise, the decoder behaves as expected in the setup, its performance decreases when the number of unreliable sensors increases. Again, in setup 2, when the number of failed sensors increases, the performance of the decoder decreases in terms of the rate of false positive and false negative. However, when a drifted sensor is detected by the decoder, we observe that the average detection delay of detection is slightly better when the number of unreliable sensors increases.

APPENDIX

Proof.

Step 1 (Necessary condition): We first prove the necessary condition for m * to be a local optimum of (3). If m < m * i , then from (4) we can decomposed

m * i n=m +1 I j=1 e ij 1 m * j ≥n (w n ij -λ) into the sum of two terms m * i n=m +1 I j∈c(i) e ij 1 m * j ≥n (w n ij -λ) and + m * i n=m +1 I j =c(i) e ij 1 m * j ≥n (w n ij - λ). From assumption A, note that for all n ≤ m * i , I j∈C(i) 1 m * j ≥n (w n ij -λ) ≤ A(i, n)(Φ 1 -λ), we have m * i n=m +1 I j∈C(i) αM ij 1 m * j ≥n (w n ij -λ) ≤ (m * i -m -1) | C(i) | αC M E (c 3 -λ), m * i n=m +1 I j =C(i) αM ij 1 m * j ≥n (w n ij -λ) ≤ (m * i -m -1) j =i | C(j) | αC M E (c 3 -λ)
Therefore it follows that if (5), then We will use a contradiction argument. We need to prove that it exists m , such that (4) is not satisfied for m loc . Two cases need to be solved. First, for a given i, let us assume that April 23, 2020 DRAFT m loc i < m ≤ m * i . Note that i ∈ I 1 (m ) and from (4) we have:

Figure 1 .

 1 Figure 1. Setup (1): The measurements which did not drift are in grey in this figure. The drift measurements are in red, purple and yellow. The vertical dotted lines represented the instant at which sensors drifted and the vertical dashed lines capture the instant predicted by the decoder.

Figure 2 .Figure 3 .

 23 Figure 2. Setup (2): The measurements which did not drift are in grey and dark blue in this figure. The measurements in dark blue indicate that the decoder predicted, wrongly, that these measurements are unreliable. The drift measurements are in green, light grey and yellow. As in figure 1, the vertical dotted lines represented the instant at which sensors drifted and the vertical dashed lines capture the instant of drift beginning predicted by the decoder.

e

  ij 1 m * j ≥n (w n ij -λ) ≤ 0. A similar line of thought can be used for every m > m * i . Instead of find an upper bound we need to find a lower bound for m n=m * i +1 I j=1 e ij 1 m * j ≥n (w n ij -λ). Step 2 (Uniqueness condition): Let us assume that there exists another local solution m loc = m * .

(( 1 -(Φ 1 e ij 1

 111 α) + αM ij )1 m loc j ≥n (w n ij -λ)αM ij 1 m loc j ≥n (w n ij -λ)From assumption A, fromB(i, n) ≤ B c(i) and from | A(i, n) |=| C(i) | -| B(i, n) | we caneasily deduce the following inequality:-λ) | A(i, n) | +(c 2 -λ) | B(i, n) | ≤ (1 -α)(m -m loc i -1) ×((Φ 1 -λ) | A(i, n * ) | +(c 2 -λ) | B(i, n * )) |,wheren * = argmax n (Φ 1 -λ) | A(i, n) | +(c 2 -λ) | B(i, n) |.Moreover, by using assumption A, we have:αM ij 1 m loc j ≥n (w n ij -λ) ≤ (m -m loc i -1) | C(i) | αC M E (c 3 -λ) αM ij 1 m loc j ≥n (w n ij -λ) ≤ (m -m loc i -1) j =i | C(j) | αC M E (c 3 -λ)Note that if the assumption of the theorem are respected then it is leading us to a contraction. m loc j ≥n (w n ij -λ) < 0 and not be positive. The second case is the following, for a given i, let us assume that m * i < m ≤ m loc i . Note that i ∈ I 2 (m ). Then by using the same argument as previously and by replacing(Φ 1 -λ) | A(i, n) | +(c 2 -λ) | B(i, n) | by (Φ 2 -λ) | A(i, n) | +(c 3 -λ) | B(i, n) |, we will reach the same conclusion. April 23, 2020 DRAFT


Pierre Fortin
email: pierre.fortin@sorbonne-universite.fr

Ambroise Fleury
email: ambroise.fleury.etu@univ-lille.fr

François Lemaire
email: francois.lemaire@univ-lille.fr

Michael Monagan
email: mmonagan@cecm.sfu.ca

High performance SIMD modular arithmetic for polynomial evaluation

Keywords: modular arithmetic, SIMD, polynomial evaluation, operational intensity, instruction-level parallelism

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Computer Algebra, also called symbolic computation, consists of developing algorithms and data structures for manipulating mathematical objects in an exact way. Multivariate polynomials with rational coefficients are essential objects in Computer Algebra, and they naturally arise in many applications (Mechanics, Biology, . . . ), especially in non-linear problems. Among the classical operations on multivariate polynomials (sum, product, quotient, . . . ), two non trivial operations are essential: polynomial factorization and polynomial gcd (greatest common divisor) computation [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Keith | Algorithms for Computer Algebra[END_REF][START_REF] Von | Modern Computer Algebra[END_REF]. Those operations are necessary for solving polynomial systems and simplifying quotients of multivariate polynomials.

Modern algorithms [START_REF] Monagan | The complexity of sparse Hensel lifting and sparse polynomial factorization[END_REF][START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF] for factorization and gcd computation rely on many polynomial evaluations, which dominate the overall computation cost. These polynomial evaluations have two main features. First, these are partial evaluations in the sense that not all variables are evaluated: given a polynomial with n variables, we evaluate n -2 variables which results in a polynomial with 2 variables. Second, the variables are evaluated at integers modulo a prime p thus all integer arithmetic operations are performed modulo p. Computing modulo a 64 bit prime p makes it possible to use machine integers and native CPU operations, instead of arbitrary-precision integers. Since these partial modular polynomial evaluations are currently a performance bottleneck for polynomial factorizations and gcd computations, we aim in this article to speed-up their computation on modern CPUs.

We focus here on one compute server since most symbolic computations are usually performed on personal workstations. We distinguish three main techniques to obtain performance gain on current CPU compute servers [START_REF] Hennessy | Computer Architecture: A Quantitative Approach, Sixth Edition[END_REF]:

• increasing the compute efficiency. This can be achieved by increasing the operational intensity, i.e. the number of operations per byte of memory (DRAM) traffic, which helps exploit modern HPC architectures whose off-chip memory bandwidth is often the most performance constraining resource [START_REF] Williams | Roofline: An insightful visual performance model for multicore architectures[END_REF]. The compute efficiency can also be increased by better filling the pipelined floating-point units with instruction-level parallelism;

• exploiting data-level parallelism on vector (or SIMD -single instruction, multiple data) units. Such parallelism is increasingly important in the overall CPU performance since the SIMD vector width has been constantly increasing from 64 bits (MMX[8] and 3DNow! [START_REF]3DNow! Technology Manual[END_REF]) to 128 bits (SSE [START_REF]Intel. Intel SSE4 Programming Reference[END_REF], AltiVec [START_REF] Diefendorf | Altivec extension to Power PC accelerates media processing[END_REF]), then to 256 bits (AVX and AVX2 [START_REF]Intel. Intel Architecture Instruction Set Extensions Programming Reference[END_REF]), and recently to 512 bits (AVX-512 [START_REF]Intel. Intel Architecture Instruction Set Extensions Programming Reference[END_REF]). For 64-bit integers, such AVX-512 SIMD units can now offer a 8x speedup with respect to a scalar computation. But efficiently exploiting such parallelism requires "regular" algorithms where the memory accesses and the computations are similar among the lanes of the SIMD vector unit. Moreover, multiple programming paradigms can be used for SIMD programming: intrinsics, compiler directives, automatic vectorization, to name a few. For key applications, it is important to determine which programming paradigm can offer the highest programming level (ease of programming and portability) without sacrificing (too much) performance.

• exploiting thread-level parallelism on the multiple cores and on the multiple processors available in a compute server. This can be achieved thanks to multi-process, multi-thread or multi-task programming.

Thread-level parallelism has already been introduced for the partial modular polynomial evaluation (see e.g. Hu and Monagan [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF]). We will therefore focus here on single-core optimizations, namely increasing the compute efficiency and SIMD parallelism, and we present in this article the following contributions.

• Multiple algorithms have been designed to efficiently compute modular multiplications in a scalar mode. We first justify why the floating-point (FP) based algorithm with FMA (fused multiply-add) of van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF] is best suited for current HPC architectures, especially regarding SIMD computing. We show that the optimized AVX version implementation of van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF] can safely be used in our polynomial evaluation, and we then propose the first (to our knowledge) implementation of such modular multiplication algorithm on AVX-512, as well as the corresponding FP-based modular addition. With respect to the reference polynomial evaluation implementation of Monagan and coworkers [START_REF] Monagan | The complexity of sparse Hensel lifting and sparse polynomial factorization[END_REF][START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF], which relies on scalar integer-based algorithms for modular operations, we detail the performance gains of our SIMD FP-based modular arithmetic for modular operation microbenchmarks and for the polynomial evaluation.

• We carefully compare intrinsics and OpenMP compiler directives for programming SIMD modular arithmetic and for their integration within the polynomial evaluation. We detail the relevance, the advantages, the issues and the performance results of both programming paradigms.

• We show how to significantly improve the performance of the modular polynomial evaluation by increasing the operational intensity via data reuse, and by filling the pipelines of the floating-point units. This is achieved thanks to the introduction of multiple "dependent" and "independent" evaluations and loop unrolling. We also show that close-to-optimal performance can be obtained without extra memory requirements.

In the rest of this article, we first introduce in Sect. 2 our partial polynomial evaluation. Then we detail in Sect. 3 our integration of SIMD computing first in modular arithmetic, then in the polynomial evaluation. Finally, we show how we can increase the compute efficiency of our polynomial evaluation in Sect. 4.

Modular polynomial evaluation in two variables

Presentation

We are given a multivariate polynomial f ∈ Z[x 1 , . . . , x n ] and we want to evaluate the variables x 3 , . . . , x n in f at integers modulo a prime p. Let F p denotes the finite field of integers modulo a prime p. The prime p is chosen so that all integer arithmetic in 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 β 3 , β 4 , . . . , β n uniformly at random from F p and compute a sequence of T evaluations

b t (x 1 , x 2 ) = f (x 1 , x 2 , β t 3 , β t 4 , . . . , β t n ) for 1 ≤ t ≤ T.
The values b t are polynomials in F p [x 1 , x 2 ]. We call them bivariate images of f . For convenience we will use the notation

β = (β 3 , β 4 , . . . , β n ) ∈ F n-2 p
, and

β t = (β t 3 , β t 4 , . . . , β t n ) ∈ F n-2 p
, so that we may write b t (x 1 , x 2 ) = f (x 1 , x 2 , β t ). Before presenting two application examples of such computation, we first emphasize that we evaluate here at powers of β, not at T different random points in F n-2 p , since this enables one to save computations. Indeed, when Zippel [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF] introduced the first sparse polynomial interpolation algorithm and used his algorithm to compute a gcd of two polynomials f and h in Z[x 1 , . . . , x n ], he used random evaluation points. To interpolate a polynomial with s terms his method solves one or more s × s linear systems modulo p. Using Gaussian elimination this does O(s 3 ) arithmetic operations in F p and requires space for O(s 2 ) elements of F p . Then Zippel [START_REF] Zippel | Interpolating polynomials from their values[END_REF] showed that if powers β 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 F p and O(s) space. In Sect. 2.3, we will see that using powers of β also reduces the polynomial evaluation cost.

Second, we also emphasize that we evaluate at β t for 1 ≤ t ≤ T , and not for 0 ≤ t < T , since the evaluation point β 0 = (1, 1, . . . , 1) may not be usable. For example, consider the following gcd problem in

Z[x 1 , x 2 , x 3 , x 4 ]. Let f = gc and h = gc + gd where c = x 1 x 2 + x 3 x 4 , d = x 3 -x 4 and g ∈ Z[x 1 , x 2 , x 3 , x 4 ]. Since gcd(c, d) = 1 we have gcd(f, h) = gcd(gc, gc + gd) = gcd(gc, gd) = g gcd(c, d) = g. But suppose we use β 0 = (1, 1). Then since d(x 1 , x 2 , 1, 1) = 0 we have gcd(f (x 1 , x 2 , β 0 ), h(x 1 , x 2 , β 0 )) = g(x 1 , x 2 , 1, 1)c(x 1 , x 2 , 1, 1) = g(x 1 , x 2 , 1, 1)(x 1 x 2 + 1).
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 β at random from F n-2 p and evaluating at β t for t starting at 1.

Application examples

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.

Polynomial factorization

Regarding polynomial factorization [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Keith | Algorithms for Computer Algebra[END_REF], Monagan and Tuncer [START_REF] Monagan | Sparse multivariate polynomial factorization: a highperformance design and implementation[END_REF][START_REF] Monagan | The complexity of sparse Hensel lifting and sparse polynomial factorization[END_REF] reduce factorization of a multivariate polynomial f in Z[x 1 , x 2 , . . . , x n ] to (i) evaluating f (x 1 , x 2 , β t ) for 1 ≤ t ≤ T , (ii) doing a computation with the bivariate images, and (iii) recovering the factors of f using sparse interpolation techniques [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF][START_REF] Zippel | Interpolating polynomials from their values[END_REF][START_REF] Ben | A deterministic algorithm for sparse multivariate polynomial interpolation[END_REF]. See Roche [START_REF] Daniel | What can (and can't) we do with sparse polynomials[END_REF] for a recent discussion on sparse polynomial interpolation methods and an extensive bibliography. If f = r i=1 f i is the factorization of f over Z then usually the factors f i have a lot fewer terms than their product f . Furthermore, because the method interpolates the coefficients of the factors f i from bivariate images in x 1 and x 2 , the largest coefficient is likely to have a lot fewer terms than the factor f i . Because of this, the evaluation step dominates the cost. The interpolation step, though more complicated, is cheaper because the coefficients of the factors f i being interpolated have far fewer terms than f which is being evaluated.

We wish to give an example to illustrate the numbers involved. We consider the factorization of determinants of symmetric Toeplitz matrices [START_REF] Monagan | The complexity of sparse Hensel lifting and sparse polynomial factorization[END_REF]. The m'th symmetric Toeplitz matrix T m is an m × m matrix in m variables x 0 , x 1 , . . . , x m-1 where T ij = x |i-j| . For example

T 4 =     x 0 x 1 x 2 x 3 x 1 x 0 x 1 x 2 x 2 x 1 x 0 x 1 x 3 x 2 x 1 x 0     .
The problem is to factor the polynomial det T m . For m = 14, det T m has s = 5, 165, 957 terms. It factors into two factors each with 34, 937 terms. The largest coefficient has u = 9, 705 terms. Thus u, the size of the interpolation problem, is 532 times smaller than s, the size of the evaluation problem.

Polynomial gcd

Given two polynomials f, h ∈ Z[x 1 , x 2 , . . . , x n ], to compute g = gcd(f, h), the parallel algorithm of Hu and Monagan [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF][START_REF] Hu | A fast parallel sparse polynomial GCD algorithm[END_REF] works by computing bivariate images of g modulo a prime p, namely,

g t (x 1 , x 2 ) = gcd(f (x 1 , x 2 , β t ), h(x 1 , x 2 , β t )) for 1 ≤ t ≤ T.
It then uses sparse interpolation techniques to interpolate g from the images g t (x 1 , x 2 ). Since g is a factor of f and h, the number of terms in g is typically much fewer than the number in f and h. Let us use the notation #f for the number of terms of a polynomial f . So for the gcd problem, typically #g max(#f, #h). Hu and Monagan [START_REF] Hu | A fast parallel sparse polynomial GCD algorithm[END_REF] present a "benchmark" problem where #f = 10 6 , #h = 10 6 , and #g = 10 4 . If one interpolates g from univariate images then the largest coefficient of g in x 1 has 1108 terms. If instead, as the authors recommend, one interpolates g from bivariate images, then the largest coefficient of g in x 1 and x 2 has only u = 122 terms. So for this problem, u is almost 10,000 times smaller than max(#f, #h), the size of the evaluation problem. Again, because of this, the authors found that the evaluations of the input polynomials f and h completely dominate the cost of polynomial gcd computations.

Thus two very central problems in Computer Algebra, namely, polynomial factorization and polynomial gcd computation are usually dominated by evaluations when there are many variables.

The matrix method

Let p be a prime and let f ∈ F p [x 1 , x 2 , . . . , x n ]. We may write f as

f = s i=1 a i x d i 1 x e i 2 M i (x 3 , . . . , x n ) ( 1
)
where a i ∈ F p are non-zero, d i and e i are non-negative integers and M i is a monomial in x 3 , . . . , x n . For β = (β 3 , β 4 , . . . ,

β n ) ∈ F n-2 p
, we want to compute T partial evaluations

b t (x 1 , x 2 ) = f (x 1 , x 2 , β t ) = s i=1 a i x d i 1 x e i 2 M i (β t ) for 1 ≤ t ≤ T. If we let m i = M i (β 3 , . . . , β n ) ∈ F p and M i (x 3 , . . . , x n ) = n k=3 x d ik k then we observe that M i (β t ) = n k=3 (β t k ) d ik = k=3 (β d ik k ) t = M i (β) t = m t i . Thus b t (x 1 , x 2 ) = f (x 1 , x 2 , β t ) = s i=1 a i x d i 1 x e i 2 m t i .
Now we can present the "matrix method" [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF] which relies on the powers of β to efficiently compute the T bivariate images. First we compute the monomial evaluation m i by evaluating M i (β 3 , . . . , β n ). To do this, let d k = deg(f, x k ) and let d = max n k=3 d k . We pre-compute tables of powers

β i k for 0 ≤ i ≤ d k for 3 ≤ k ≤ n.
This takes at most (n -2)(d -1) multiplications. Then, for 1 ≤ i ≤ s we compute m i = M i (β), 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

b 1 (x 1 , x 2 ) = f (x 1 , x 2 , β) = s i=1 a i m i x d i 1 x e i 2
needs 1 multiplication for each a i m i ∈ F p , hence s multiplications in total. We can compute the next evaluation

b 2 (x 1 , x 2 ) = f (x 1 , x 2 , β 2 ) = s i=1 a i m 2 i x d i 1 x e i 2 Algorithm
for i ≤ j < i + J do 7: a[j] ← a[j] ⊗ p m[j]
Hadamard product

8: c ← c ⊕ p a[j] coefficient reduction 9:
end for 10:

if c = 0 then add cx d i 1 x e i 2 to bivariate image b t 11: i ← i + J 12:
end while 13: end for using another s multiplications if we save a i m i ∈ F p for 1 ≤ i ≤ 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 t i for 1 ≤ t ≤ T, 1 ≤ i ≤ s, hence O(nd + ns + sT ) multiplications in total. With random points instead of powers of β, the T evaluations would have required a larger operation count of O(nT (d + s)) multiplications in total [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF].

One way to see the matrix method is to think of evaluating f at β t as the following t × s matrix-vector multiplication.

     m 1 m 2 . . . m s m 2 1 m 2 2 . . . m 2 s . . . . . . . . . . . . m T 1 m T 2 . . . m T s           a 1 x d 1 1 x e 1 2 a 2 x d 2 1 x e 2 2
. . .

a s x ds 1 x es 2      =      b 1 b 2 . . . b T     
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 , . . . , a s ], m = [m 1 , m 2 , . . . , m s ] and

X = [x d 1 1 x e 1 2 , x d 2 1 x e 2 2 , . . . , x ds 1 x es 2 ]. Let u•v denote the Hadamard product of two vectors u, v ∈ F s p , that is u • v = [u 1 v 1 , u 2 v 2 , . . . , u s v s ] ∈ F s p .
Then, viewing b 1 (x 1 , x 2 ) and b 2 (x 1 , x 2 ) as vectors of terms, we have

b 1 (x 1 , x 2 ) = (a • m) • X and b 2 (x 1 , x 2 ) = ((a • m) • m) • X.
Finally, for a given t we will have to compute the sum c i,t of a i m t i for all a i m t i x d i 1 x e i 2 sharing the same d i and e i values. This sum c i,t is indeed the coefficient of

x d i 1 x e i 2 in b t (x 1 , x 2
). These "coefficient reductions" are required since in Eq. ( 1), multiple M i (x 3 , . . . , x n ) can potentially share the same d i and e i values. If the monomials x d i 1 x e i 2 M i (x 3 , . . . , x n ) in the input polynomial f are sorted in lexicographical order with

x 1 > x 2 > x 3 > • • • > x n then the monomials x d i
1 x e i 2 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 , β t ) for 1 ≤ t ≤ 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 t i values in the a vector, and where ⊕ p and ⊗ p denote the ). 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 × 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.

Multi-core parallel evaluation

Hu and Monagan [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF] parallelized the matrix method for partial modular polynomial evaluations on a multi-core architecture with N cores by doing N evaluations at a time. They

first compute Γ = [m N 1 , m N 2 , . . . , m N s ] using exponentiations by squaring, requiring O(s log 2 N ) multiplications. Then they compute Λ k = a • [m k 1 , m k 2 , . . . , m k s ] for 1 ≤ k ≤ N using N s multiplications. Then, in parallel, the k'th core successively computes f (x 1 , x 2 , β k+N ) = Λ k • Γ • X; f (x 1 , x 2 , β k+2N ) = Λ k • Γ • Γ • X; f (x 1 , x 2 , β k+3N ) = Λ k • Γ • Γ • Γ • X;
. . . using Algorithm 1 (with each Λ k as the a vector, and Γ as the m vector). This was implemented with multi-task programming in Cilk [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF]. This method significantly increases the space needed as N vectors C 1 , C 2 , . . . , C N of length s are required where s can be very large. Monagan and Tuncer [START_REF] Monagan | Sparse multivariate polynomial factorization: a highperformance design and implementation[END_REF] introduced an alternative parallelization strategy by using a 1D block decomposition of the a and m vectors for each evaluation.

Finally, we mention the asymptotically fast method [START_REF] Van | On the Bit-complexity of Sparse Polynomial and Series Multiplication[END_REF] for computing f (x 1 , x 2 , β t ). If p is chosen of the form p -1 = 2 k q with 2 k > T so that an FFT of order 2 k can be done in the finite field F p , after computing the monomial evaluations Monagan and Wong [START_REF] Monagan | Fast Parallel Multi-point Evaluation of Sparse Polynomials[END_REF] found that a serial implementation of this fast method first beat the matrix method at T = 504 but that it was much more difficult to parallelize than the matrix method -the fast method required s 10 6 to deliver good parallel speedups. Moreover, the simplicity and data locality of the matrix method makes it very suitable for vectorization and other single-core optimizations targeted in this article.

m = [m 1 , m 2 , . . . , m s ], this method computes f (x 1 , x 2 , β t ) for 1 ≤ t ≤ T in O(s log 2 T ) multiplications.

SIMD modular arithmetic for partial polynomial evaluation 3.1 Selection of the modular arithmetic algorithm

Given a fixed1 integer p > 1, we focus on the efficient computation of c = a⊗ p b and d = a⊕ 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 ⊕ p can be implemented with a compare instruction, ⊗ p requires integer divisions which are expensive operations on current processors [START_REF] Fog | Optimizing software in c++. an optimization guide for windows, linux and mac platforms[END_REF], 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 ⊗ p . We briefly recall the most important ones.

In order to compute c = a⊗ p b, one can first rely on floating-point arithmetic to compute q = a×b p and then deduce c = a × b -q × p (see for example Alverson [START_REF] Alverson | Integer division using reciprocals[END_REF], Baker [START_REF] Henry | Computing a*b (mod n) efficiently in ansi c[END_REF]). This requires conversions to/from floating-point numbers, and the number of bits of the floatingpoint number mantissa has to be twice as large as the number of bits of p (to hold the product).

Algorithm 2 Modular multiplication of 64-bit integers x and y with a 50-bit prime p. x and y are considered to be already reduced modulo p, and converted to double along with p prior to the beginning of the algorithm. u stores: 1/(double) p

1: double h ← x * y ; 2: double ← fma(x, y, -h) ; 3: double b ← h * u ; 4: double c ← floor(b) ; c is the quotient ±1 5: double d ← fma(-c, p, h) ; 6: double g ← d + ;
g is the remainder ±p 7: if g ≥ p then return g -p ; 8: if g < 0.0 then return g + p ; 9: return g ; Algorithm 3 Modular addition of 64-bit integers x and y with a 50-bit prime p. x and y are considered to be already reduced modulo p, and converted to double along with p prior to the beginning of the algorithm.

1: double s ← x + y ; 2: return s ≥ p ? s -p : s ;

In order to avoid the conversions between floating point numbers and integers, the floatingpoint 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 integerbased method is known as the Barrett's product [START_REF] Barrett | Implementing the rivest shamir and adleman public key encryption algorithm on a standard digital signal processor[END_REF][START_REF] Granlund | Division by invariant integers using multiplication[END_REF]. Another integer-based approach relies on Montgomery's reduction [START_REF] Peter | Modular multiplication without trial division[END_REF]: a comparison between the two methods has been done for example by van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF]. An improved version of Barrett's product with integer only operations has been proposed by Möller and Granlund [START_REF] Moller | Improved division by invariant integers[END_REF]: Monagan and coworkers [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF][START_REF] Monagan | Sparse multivariate polynomial factorization: a highperformance design and implementation[END_REF] use an implementation (written by Roman Pearce) of this latter [START_REF] Moller | Improved division by invariant integers[END_REF] method (with p -1 precomputed) in their original code for polynomial evaluation with 64-bit integers. This offers a 11x performance gain [START_REF] Gibson | Optimizing and Parallelizing the Modular GCD Algorithm[END_REF] for ⊗ 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 [START_REF] Moller | Improved division by invariant integers[END_REF] 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 [START_REF] Ogita | Accurate sum and dot product[END_REF]. Such error-free transformation computes the accurate floating-point result of the product. As described and proved by van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF], 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 [START_REF] Ogita | Accurate sum and dot product[END_REF], and hence to precisely handle the multiplication result before reduction modulo p. More precisely, stores the rounding error of the product x * y (i.e. h + 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 in g in order to take into account the initial rounding error. g is finally corrected so that we exactly have: g ≡ x × y (mod 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 , . . . , x n ) = s i=1 a i M i (x 1 , x 2 , . . . , x n ) be a polynomial we wish to interpolate. Ben-Or/Tiwari [START_REF] Ben | A deterministic algorithm for sparse multivariate polynomial interpolation[END_REF] and Zippel [START_REF] Zippel | Interpolating polynomials from their values[END_REF] both pick β ∈ F n p at random and interpolate f from f (β), f (β2 ), f (β 3 ), . . . . Both methods require the monomial evaluations to be distinct, that is, M i (β) = M j (β) for 1 ≤ i < j ≤ 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.

SIMD programming paradigms

We plan to integrate the SIMD implementations of the FMA-based ⊕ p and ⊗ 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 2 (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.

SIMD intrinsics and the AVX-512 version

Using AVX intrinsics: Van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF] have presented a SSE/AVX version of Algorithm 2 to implement ⊗ p . They use two SSE/AVX blendv pd intrinsics to efficiently implement the two final tests (Lines 7-8 in Algorithm 2), hence removing divergence in the SIMD computation. This blendv pd intrinsic blend double elements from two vectors depending on the most significant bit of elements from a third vector. For floating-point elements this most significant bit corresponds to the sign bit, which enables one to implement in SIMD without branching the two final tests using comparisons to 0.0. We will rely on this AVX version on AVX2 CPUs. However, we recall that IEEE standard 754 for floating-point arithmetic [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF] includes signed zeros which may lead to incorrect results regarding the use of blendv pd. Indeed if for example g equals -0.0 at Line 8 in Algorithm 2, using the blendv pd instruction directly on g would return p, which is incorrect for modulo p arithmetic. We show below that -0.0 cannot appear in our specific context: we can thus safely use this implementation with AVX intrinsics.

Van der Hoeven et al. also present an AVX ⊕ p implementation (see function 3.9 [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF]) of the scalar FP-based ⊕ p algorithm (presented in Algorithm 3). Similarly, as shown below, we can safely ignore signed zeros for the the blendv pd instruction used in this AVX version.

Regarding the issue with signed zeros and the AVX blendv pd intrinsic: For ⊗ p , we consider the tests at Lines 7 (g ≥ p here rewritten as g -p ≥ 0) and 8 (g > 0) in Algorithm 2. We aim at showing here that no -0.0 value will occur in these two tests, so that one can safely use the blendv pd intrinsics to implement in AVX the conditional affectation resulting from these tests. For completeness, we recall beforehand the paragraphs 3 and 4 of §6.3, The sign bit, from the IEEE 754 standard [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Std[END_REF].

Paragraph 3: When the sum of two operands with opposite signs (or the difference of two operands with like signs) is exactly zero, the sign of that sum (or difference) shall be +0 in all rounding-direction attributes except roundTowardNegative; under that attribute, the sign of an exact zero sum (or difference) shall be -0. However, x + x = x -(-x) retains the same sign as x even when x is zero.

Paragraph 4: When (a * b) + c is exactly zero, the sign of fusedMultiplyAdd(a, b, c) shall be determined by the rules above for a sum of operands. When the exact result of (a * b) + c is non-zero yet the result of fusedMultiplyAdd is zero because of rounding, the zero result takes the sign of the exact result. Now, let us first consider the case where x or y is zero (possibly both). Then h computed at Line 1 is either +0.0 or -0.0. According to the paragraph 4 quoted just above, the computed at Line 2 is obtained from the computation h -h which cannot result in -0.0 thanks to paragraph 3 and since we rely on the default rounding mode (Round to nearest). Thus = +0.0. Using again paragraph 3, the g computed at Line 6 by g = d + cannot be equal to -0.0 since = +0.0. Therefore g at Line 6 equals +0.0 (since the expected result of the modular product is zero here): thus no -0.0 is evaluated in the two tests at Lines 7 and 8.

Second, consider the case were both x and y are nonzero. Since p is a prime number, and since x < p and y < p, the product x × y cannot be zero modulo p: this is due to the uniqueness of the prime factorization of x × y. Thus the function returns a nonzero value that is not a multiple of p. Consequently g cannot hold -0.0, and neither can g -p evaluate to -0.0.

For ⊕ p , the blendv pd intrinsic evaluates the result of a subtraction by p since the test s ≥ p at Line 2 in Algorithm 3 is rewritten as s -p ≥ 0 (see function 3.9 [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF]). According to paragraph 3, this subtraction cannot lead to -0.0 since p = -0.0. Algorithm 4 AVX-512 modular multiplication of AVX-512 vectors x and ȳ with a 50-bit prime p replicated in the AVX-512 vector p. Elements of x and ȳ are considered to be already reduced modulo p, and converted (like p) to double elements (x, ȳ, p and ū being m512d vectors) prior to the beginning of the algorithm. ū stores replicates of: 1/(double) p 1: m512d h ← mm512 mul pd(x, ȳ) ; 2: m512d ¯ ← mm512 fmsub pd(x, ȳ, h) ; 3: m512d b ← mm512 mul pd( h, ū) ; 4: m512d c ← mm512 floor pd( b) ; 5: m512d d ← mm512 fnmadd pd(c, p, h) ; 6: m512d ḡ ← mm512 add pd( d, ¯ ) ; 7: mmask8 m ← mm512 cmplt pd mask(ḡ, mm512 setzero pd()) ; 8: mmask8 mm ← mm512 cmple pd mask(p, ḡ) ; 9: ḡ ← mm512 mask add pd(ḡ, m, ḡ, p) ; 10: ḡ ← mm512 mask sub pd(ḡ, mm, ḡ, p) ; 11: return ḡ ; Algorithm 5 AVX-512 modular addition of AVX-512 vectors x and ȳ with a 50-bit prime p replicated in the AVX-512 vector p. Elements of x and ȳ are considered to be already reduced modulo p, and converted (like p) to double elements (x, ȳ and p being m512d vectors) prior to the beginning of the algorithm.

1: m512d s ← mm512 add pd(x, ȳ) ; 2: mmask8 m ← mm512 cmple pd mask(p, s) ; 3: return mm512 mask sub pd(s, m, s, p) ;

Using AVX-512 intrinsics: Regarding AVX-512, the blendv pd intrinsic is not available: we therefore explicitly build 8-bit masks to conditionally perform (without branching) the addition and the subtraction at the end of the algorithm, as presented in Algorithm 4. Hence the SIMD divergence is efficiently handled within the AVX-512 arithmetic instructions. To our knowledge this is the first AVX-512 floating-point based modular arithmetic. Orisaka et al. [START_REF] Orisaka | XVIII Simpósio Brasileiro de Segurança da Informação e Sistemas Computacionais[END_REF] have also accelerated modular arithmetic with AVX-512 but using Montgomery reduction and targeting very large primes for cryptography.

Regarding ⊕ p , we also adapt the SSE/AVX implementation for floating-point numbers presented by Van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF] to AVX-512. As detailed in Algorithm 5, we use 1 addition, 1 comparison and 1 masked subtraction in AVX-512 instead of 1 addition, 1 subtraction and 1 blendv pd in AVX.

Microbenchmarks

We start with microbenchmarks, presented in Fig. 1, of the two modular arithmetic operations ⊗ p and ⊕ 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 [START_REF] Gibson | Optimizing and Parallelizing the Modular GCD Algorithm[END_REF] on older CPUs.

Regarding ⊗ 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 implementa- tion 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 ⊗ p implementations is indeed twice greater than the AVX2 one. Regarding ⊕ 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 integerbased ⊕ p implementation can indeed be replaced by shifting [START_REF] Law | A parallel implementation for polynomial multiplication modulo a prime[END_REF][START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF], 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 ⊕ p implementation is twice faster than the AVX2 one of Van der Hoeven et al. [START_REF] Van Der Hoeven | Modular SIMD arithmetic in Mathemagix[END_REF]. 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 3 , 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 ⊗ 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 ⊗ 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 ⊗ 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 ≥ 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 ⊕ 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 ⊕ 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).

Integration in polynomial evaluation

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 Algorithm 6 SIMD compute kernel of the matrix method (see Algorithm 1 for notations and inputs). x denotes the SIMD vector corresponding to variable x. V is the size of the SIMD vector.

1: In-place conversions for vectors a and m (64-bit integers → doubles) 2: for each evaluation 1 ≤ t ≤ T do

3: i ← 1; b t ← 0 4: while i ≤ s do 5: c ← 0.0 6:
J ← #monomials with same (d i , e i )

7:

for i ≤ j < i + J with step V do 8: ā ← a[j .. j + V -1] SIMD load 9: m ← m[j .. j + V -1] SIMD load 10: ā ← ā ⊗ p m SIMD Hadamard product 11: c ← c ⊕ p ā SIMD coefficient reduction 12: a[j .. j + V -1] ← ā SIMD store 13:
end for if c = 0.0 then convert c to 64-bit integer and add cx d i 1 x e i 2 to bivariate image b t 16:

i ← i + J 17:
end while 18: end for 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 [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF] (using integer-based modular arithmetic). This is due to the slow FP-based ⊕ 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 ⊕ p and its branch instruction, we use SIMD shuffle instructions to write a parallel tree-shaped reduction using only SIMD ⊕ 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 × 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 ⊕ 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 ⊗ p and ⊕ 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 ⊗ p (see Algorithm 2) and 2 flop for ⊕ 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.

Increasing the compute efficiency

We now focus on the compute efficiency of our SIMD polynomial evaluation. More precisely, we aim to fill at best the pipelined floating-point units and to minimize the time lost in memory accesses.

Multiple dependent evaluations

We first rely on the consecutive powers of β 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 ⊕ p operation of the t evaluation and the ⊗ 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.

Multiple independent evaluations

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 [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF] 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

Γ = [m T i 1 , m T i 2 , . . . , m T i s ] using O(s log 2 T i ) SIMD multiplications. We also precompute Λ k = a • [m k+1
1 , m k+1 2 , . . . , m k+1 s ] for 0 ≤ 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 (Λ k ) 0≤k<T i ), then the V is the size of the SIMD vector. For the ease of reading we only represent here the Hadamard products, but the SIMD coefficient reductions are performed alongside, and the ck i ,k d final reductions when required. First (Fig. 3(a)), we precompute Γ and (Λ k ) 0≤k<T i =3 . Then (Fig. 3(b)), using Γ we compute in SIMD the first V products of the first

Λ 0 = a 1 m 1 . . . a i m i . . . a s m s Λ 1 = a 1 m 2 1 . . . a i m 2 i . . . a s m 2 s Λ 2 = a 1 m 3 1 . . . a i m 3 i . . . a s m 3 s Γ = m 3 1 . . . m 3 i . . . m 3 s (a) Precomputations Λ 0 = a 1 m 4 1 . . . a V m 4 V . . . a s m 4 s Λ 1 = a 1 m 5 1 . . . a V m 5 V . . . a s m 5 s Λ 2 = a 1 m 6 1 . . . a V m 6 V . . . a s m 6 s Λ 0 = a 1 m 7 1 . . . a V m 7 V . . . a s m 7 s Λ 1 = a 1 m 8 1 . . . a V m 8 V . . . a s m 8 s Λ 2 =
T i × T d = 3 × 2 evaluations in (Λ k ) 0≤k<3 .
The next V products are performed in SIMD for the same evaluations (Fig. 3(c)), and so on for the remainings of vectors (Λ k ) 0≤k<3 and Γ. Once these first evaluations have been fully computed, we start again with the next 3 × 2 evaluations (Fig. 3(d)), until all evaluations have been fully processed.

Hadamard product Λ k ← Λ k • Γ (with 0 ≤ 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 ⊕ p and ⊗ 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 Λ k of the complete a vector. Moreover, for each independent evaluation k we have to load the Λ k vector from memory and store its update in memory. The operational intensity is thus only improved for the Γ 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 × 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 × T d , the remaining evaluations are processed first by blocks of T i × N d evaluations (with 1 ≤ N d < T d ) and then with N d dependent evaluations (with N d < T i ), where T mod

(T i × T d ) = T i × N d + N d .
The final version is presented in Algorithm 7, along with the SIMD programming. Once Algorithm 7 SIMD compute kernel of the matrix method with independent and dependent evaluations (see Algorithm 1 for notations and inputs). x denotes the SIMD vector corresponding to variable x, and xi the ith SIMD vector. V is the size of the SIMD vector.

1: Pre-compute Γ = [m T i 1 , . . . , m T i s ] and Λ k i = a • [m k i +1 1 , . . . , m k i +1
s ] for 0 ≤ k i < T i (with in-place conversions: 64-bit integers → doubles) 2: for each evaluation 1 ≤ t ≤ T with step T i × T d do

3: i ← 1; b t+k i +k d T i ← 0, for (0 ≤ k i < T i ; 0 ≤ k d < T d ) 4: while i ≤ s do 5: ck i ,k d ← 0.0, for (0 ≤ k i < T i ; 0 ≤ k d < T d ) 6:
J ← #monomials with same (d i , e i )

7:

for i ≤ j < i + J with step V do 8:

Λk i ← Λ k i [j .. j + V -1], for 0 ≤ k i < T i SIMD loads 9: Γ ← Γ[j .. j + V -1] SIMD load 10: for 0 ≤ k d < T d do 11:
for 0 ≤ k i < T i do

Λ k i [j .. j + V -1] ← Λk i , for 0 ≤ k i < T i SIMD stores 17:
end for 18:

for 0 ≤ k d < T d do 19:
for 0 ≤ k i < T i do 20: i ← i + J

c k i ,k d ← reduce(c k i ,k d , ⊕ p ) ck i ,k d final reduction 21: if c k i ,k d = 0.0 then convert c k i ,

25:

end while 26: end for T i × T d evaluations have been computed together for some monomials, we could choose to iterate over the next T i × T d evaluations or to iterate over the next monomials with same (d i , e i ) values. If we had iterated over the next T i × T d evaluations, we would have had to store

T i × T d SIMD vectors (all (c k i ,k d ) 0≤k i <T i , 0≤k d <T d )
for the coefficient reductions. By iterating on the next monomials (as done in Algorithm 7), T i + 1 SIMD loads (all ( Λk i ) 0≤k i <T i and for Γ) and T i SIMD stores (all ( Λk i ) 0≤k i <T i ) are required. As T d > 2 in practice (as confirmed for the optimal configurations in Sect. 4.4), it is indeed preferable to iterate over the next monomials in order to minimize the number of memory accesses and hence increase the operational intensity. This results in the end in an algorithm where we browse all the monomials with same (d i , e i ) values to compute T i × T d evaluations at a time. An illustration of the execution of Algorithm 7 is given in Fig. 3.

Loop unrolling

Loop unrolling [START_REF] Hennessy | Computer Architecture: A Quantitative Approach, Sixth Edition[END_REF] enables us to remove the exit test at the end of the loop body and to interleave instructions from successive loop iterations in order to better fill the pipelines. The two nested loops over the T d and T i evaluations (Lines 10 and 11 in Algorithm 7) are hence completely unrolled thanks to the "unroll(F)" pragma of ICC and to the "GCC unroll F" pragma of GCC (recently introduced in GCC 8) to impose an unroll factor of F (F being respectively equal to T d and T i ). Similarly, we unroll the third loop over the monomials with same (d i , e i ) values (Line 7 in Algorithm 7) by a factor M . We could also have let the compiler choose which loops to unroll (or not) and determine the best unroll factors. This leads to similar performance with GCC, but to lower performance with ICC (up to 8.5% performance loss). We thus impose our unrollings on the three loops with the corresponding pragmas.

Once the three nested loops have been unrolled, we rely on the compiler and on the out-oforder execution of the processor to schedule at best the instructions to fill the pipelines and to overlap the memory accesses. Other loops over T d and/or T i evaluations (Lines 8, 16, 18, 19, in Algorithm 7) are also similarly unrolled.

Performance results

Using the test platforms (server and compiler) described in Table 1, we present in Figure 4 the performance results for all possible configurations for (T i , T d , M ), each value ranging in 1,2,4,8,16. We also indicate the performance of the Base SIMD code corresponding to the version obtained with the SIMD intrinsics only (as presented in Sect. 3). The performance varies significantly depending on the (T i , T d , M ) values, especially on T i and T d , which shows the relevance of these parameters. The performance impact of M is lower but can still reach 11% for some (T i , T d ) configurations. The best configurations are the following.

• (T i = 8, T d = 16, M = 1) with GCC on AVX2 server: 15.65 Gflop/s, and 42% of As determining the theoretical peak performance of modern CPUs becomes more and more complicated [START_REF] Dolbeau | Theoretical peak FLOPS per instruction set: a tutorial[END_REF], we use the BLAS DGEMM routine of the Intel MKL 4 to estimate the singlecore double-precision peak performance at 45 Gflop/s on AVX2 server and 101 Gflop/s on AVX-512 server. Moreover, we can only reach 61% of the peak performance since there are only 2 FMAs out of the 9 floating-point instructions required for ⊗ p and ⊕ p . We manage hence to reach 57% and 63% of the attainable single-core peak performance, respectively on AVX2 server and on AVX-512 server.

In the end, as shown in Table 2, we manage to reach speedups of almost 5x and 10x (respectively on AVX2 server and on AVX-512 server) on one CPU core over the reference original polynomial evaluation with scalar integer-based modular arithmetic. It can be noticed that Monagan and coworkers already used to process two (dependent) evaluations at a time to increase the operational intensity for some variants of the polynomial evaluation. However for the variant studied in this article (see Algorithm 1), which is the fastest one, no performance gain is obtained by processing two evaluations at a time with the original scalar code. Such divergence with respect to the gains obtained in Fig. 4 can be explained by the lack of other optimizations (multiple independent evaluations, loop unrollings) as well as by the differences in the modular arithmetic implementation between the original integer-based version by Monagan and coworkers and the floating-point based version of this article.

Without extra memory requirements

One drawback of using multiple independent evaluations is the significant memory overhead: the complete a vector has to be duplicated for each extra independent evaluation. We therefore investigate here the best attainable performance without any extra independent evalution.

We first implement a code without multiple independent evaluation, and tune the T d and M parameters for this code via extensive benchmarks (as in Sect. 4.4). Figure 5 shows for each test platform the performance drop obtained for this code (referred to as Without T i ) with respect to the best version obtained in Sect. 4.4 (referred to as With T i ). The performance drop is important here (up to 36%), due to the lower number of independent instructions to fill the pipelines.

We then introduced manual loop unrolling, using preprocessor macros to ease and automate the tedious code writing. We also manually group all arithmetic instructions. This way, we provide all arithmetic instructions for the computation of T d dependent evaluations and M × V monomials to the compiler and to the out-of-order execution of the processor, so that these can be scheduled at best to fill the pipelines. One can see that this new version (referred to as: Without T i , manual unrolling; and after tuning of its T d and M parameters) greatly reduces the performance drop with respect to the best version obtained in Sect. 4.4 (With T i ). This way, we can reach 95% (GCC) and 94% (ICC) on AVX2 server and 93% (GCC) and 98% (ICC) on AVX-512 server of the best attainable performance (With T i ). At the price of non-negligeable development efforts, we can thus obtain, without introducing extra memory, almost the same performance of our best versions with multiple independent evaluations. It can also be noticed that the performance impact of M is here much more important than in Sect. 4.4 (detailed tests not shown). Such manual loop unrolling and instruction grouping have also been tried on the best version obtained in Sect. 4.4 (with multiple independent evaluations): this however only offers up to 5.1% performance gain for such code. In our opinion, this does not justify the manual unrolling development effort when using multiple independent evaluations.

Conclusion

In this article, we have first justified the choice of a modular multiplication algorithm relevant for HPC and SIMD computing. We have ensured the correct use of an optimized AVX2 implementation (regarding a potential issue with signed zeros and the blendv pd intrinsic) and we have presented its AVX-512 version. This floating-point (FP) based algorithm with FMAs (fused multiply-adds) enables us to obtain SIMD speedups of up to 3.7x on AVX2, and up to 7.2x on AVX-512, which validates its efficiency. With respect to a reference (scalar) integer-based modular arithmetic, the performance gains are similar for our SIMD FP-based modular multiplication and for the corresponding SIMD FP-based modular addition. As all current desktop and HPC processors have SIMD units, we believe that such SIMD FP-based modular arithmetic should from now on be used instead of the scalar ones. Using OpenMP for their SIMD programming turned out to be a very promising approach on the new AVX-512 units (with GCC), due to its very relevant performance-programmability trade-off. Currently, we still rely on intrinsics programming for best performance and performance portability among compilers.

In a second part, we have focused on the partial polynomial evaluation which is a key computation in Computer Algebra. We have rewritten this algorithm in order to introduce multiple independent and dependent evaluations. These enable us, along with loop unrolling, to fill at best the pipelined floating-point units of the CPU and to minimize the time lost in memory accesses. Combined with SIMD computing, we achieve speedups up to almost 5x on AVX2 and up to almost 10x on AVX-512 with respect to the reference implementation of the polynomial evaluation. Moreover, using manual loop unrolling we manage to closely reach such performance gains without extra memory requirements.

In the future, we plan to integrate our efficient polynomial evaluation on one CPU core in the multi-core parallel implementation of Monagan and coworkers [START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF][START_REF] Monagan | Sparse multivariate polynomial factorization: a highperformance design and implementation[END_REF], and to study the performance impact on polynomial factorizations and polynomial greatest common divisor computations. We also believe that GPUs may be well suited to further accelerate our polynomial evaluation thanks to their higher compute power and memory bandwidth. We emphasize that our FP-based modular arithmetic will be very relevant for the GPU FMA SIMD units, and will offer a direct and efficient implementation of modular arithmetic on GPUs. We may also investigate using a few less bits for our prime p in order to decrease the number of reductions as done for example with error-free transformations in linear algebra [START_REF] Jean | A parallel algorithm for dot product over word-size finite field using floating-point arithmetic[END_REF].

  arithmetic operators modulo p (c = a⊗ p b denoting c ≡ a × b (mod p), and d = a⊕ p b denoting d ≡ a + b (mod p) with (a, b, c, d) ∈ F 4 p

Figure 1 :

 1 Figure 1: Performance comparison of various implementations for ⊗ p and ⊕ p . Integer denotes the original integer-based implementation as used by Monagan and coworkers for ⊗ p[START_REF] Hu | A Fast Parallel Sparse Polynomial GCD Algorithm[END_REF][START_REF] Monagan | Sparse multivariate polynomial factorization: a highperformance design and implementation[END_REF] and ⊕ p[START_REF] Law | A parallel implementation for polynomial multiplication modulo a prime[END_REF]. FP denotes the scalar floating-point based implementation without considering the conversions between integers and floating-points numbers, whereas FP conv. includes these conversions. Intr. FP (respectively OpenMP FP) denotes the SIMD version of the floatingpoint based implementation using intrinsics (resp. using OpenMP). These performance results have been obtained on element-wise operations over vectors of 2048 elements.

14 :c

 14 ← reduce(c, ⊕ p ) c final reduction 15:

  On AVX-512 server.

Figure 2 :

 2 Figure 2: Performance comparison of scalar and SIMD polynomial evaluations. Integer denotes the original polynomial evaluation by Monagan and coworkers[5] with scalar integer-based modular arithmetic. FP denotes the use of scalar floating-point modular arithmetic for the polynomial evaluation. Intr. FP (respectively OpenMP FP) denotes the SIMD version of the floating-point based modular arithmetic using intrinsics (resp. using OpenMP).

a 1 m 9 1 . . . a V m 9 V. . . a s m 9 s= a 1 m 4 1 . . . a V+1 m 4 V+1. . . a 2V m 4 2V. . . a s m 4 sΛ 1 = a 1 m 5 1 . . . a V+1 m 5 V+1. . . a 2V m 5 2V. . . a s m 5 sΛ 2 = a 1 m 6 1 . . . a V+1 m 6 V+1. . . a 2V m 6 2V. . . a s m 6 sΛ 0 = a 1 m 7 1 . . . a V+1 m 7 V+1. . . a 2V m 7 2V. . . a s m 7 sΛ 1 = a 1 m 8 1 . . . a V+1 m 8 V+1. . . a 2V m 8 2V. . . a s m 8 sΛ 2 = a 1 m 9 1 . . . a V+1 m 9 V+1. . . a 2V m 9 2V. . . a s m 9 s= a 1 m 10 1 . . . a V m 10 V. . . a s m 10 sΛ 1 = a 1 m 11 1 . . . a V m 11 V. . . a s m 11 sΛ 2 = a 1 m 12 1 . . . a V m 12 V. . . a s m 12 sΛ 0 = a 1 m 13 1 . . . a V m 13 V. . . a s m 13 sΛ 1 = a 1 m 14 1 . . . a V m 14 V. . . a s m 14 sΛ 2 = a 1 m 15 1 . . . a V m 15 V. . . a s m 15 sFigure 3 :

 19914441155521666017771188821999110101111112112120113131114142115153 Figure3: Illustration of the execution of Algorithm 7 with T i = 3, T d = 2 and using notations of Sect. 2.3. V is the size of the SIMD vector. For the ease of reading we only represent here the Hadamard products, but the SIMD coefficient reductions are performed alongside, and the ck i ,k d final reductions when required. First (Fig.3(a)), we precompute Γ and (Λ k ) 0≤k<T i =3 . Then (Fig.3(b)), using Γ we compute in SIMD the first V products of the first T i × T d = 3 × 2 evaluations in (Λ k ) 0≤k<3 . The next V products are performed in SIMD for the same evaluations (Fig.3(c)), and so on for the remainings of vectors (Λ k ) 0≤k<3 and Γ. Once these first evaluations have been fully computed, we start again with the next 3 × 2 evaluations (Fig.3(d)), until all evaluations have been fully processed.

  With ICC on AVX-512 server.

Figure 4 :

 4 Figure 4: Performance results for all possible (T i , T d , M ) configurations.

Figure 5 :

 5 Figure5: Performance comparison between best configuration with extra memory for multiple independent evaluations (with T i ) and best configurations without extra memory for multiple independent evaluations (without T i ), with and without manual loop unrolling.

1

  Compute kernel of the matrix method for T bivariate images of a polynomial f ∈ F p [x 1 , x 2 , . . . , x n ], using notations of Eq. (1). Inputs are the vector m = [M 1 (β), . . . , M s (β)] ∈ F s p of monomial evaluations and the coefficient vector a = [a 1 , . . . , a s ] ∈ F s p . 1: for each evaluation 1 ≤ t ≤ T do

	2:	i ← 1; b t ← 0
	3:	while i ≤ s do
	4:	c ← 0
	5:	J ← #monomials with same (d i , e i )
	6:

Table 1 :

 1 Test platforms.

	Servers	Name:	Hardware features:
		AVX2 server	2× Intel Xeon CPU E5-2695 v4 CPUs: 2×18 2-way SMT
			cores -2.10 GHz (base) / 3.30 GHz (turbo) -AVX2
		AVX-512 server 2× Intel Xeon Gold 6152 CPUs:	2×22 2-way SMT
			cores -2.10 GHz (base) / 3.70 GHz (turbo) -AVX512
	Compilers Name:	Performance-related options:
		GCC 8.2.0	-O3 -mfma -fno-trapping-math -march=native -mtune=native
		ICC 19.0.3.199	-O3 -fma -xhost

12 :

 12 ck i ,k d ← ck i ,k d ⊕ p Λk i SIMD coefficient reduction

	13:	Λk i ← Λk i ⊗ p Γ	SIMD Hadamard product
	14:	end for
	15:	end for
	16:

  k d to 64-bit integer and add c k i ,k d x d i 1 x e i 2 to bivariate image b t+k i +k d T i

	22:	end for
	23:	end for
	24:

Table 2 :

 2 Performance comparison between our best version and the reference implementation using integer-based arithmetic. (T i = 8, T d = 8, M = 8) with ICC on AVX2 server: 15.21 Gflop/s, and 29% of performance gain over the Base SIMD code.• (T i = 8, T d = 8, M = 4) with GCC on AVX-512 server: 39.11 Gflop/s, and 121% of performance gain over the Base SIMD code.• (T i = 8, T d = 8, M = 1) with ICC on AVX-512 server: 37.31 Gflop/s, and 109% of performance gain over the Base SIMD code.

			Reference scalar	SIMD FP-based version with
	Server	Compiler	integer-based version	improved compute efficiency	Gain
			(time in ms)	(time in ms)
	AVX2 server	GCC	15262	3482	4.4x
	AVX2 server	ICC	17704	3671	4.8x
	AVX-512 server	GCC	12984	1411	9.2x
	AVX-512 server	ICC	14638	1476	9.9x
	performance gain over the Base SIMD code.

•

The 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.

https://www.openmp.org/

More 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.

See: https://software.intel.com/en-us/mkl

Acknowledgments

The authors would like to thank the master in computer science at Sorbonne Université, especially N. Picot and P. Cadinot, for administering and providing access to the compute servers. They also thank Professor S. Graillat (Sorbonne Université) for helpful discussions on error-free transformations.