Dominic Bucerzan

Pierre-Louis Cayrel
email: pierre.louis.cayrel@univ-st-etienne.fr

Vlad Dragoi
email: vlad.dragoi1@univ-rouen.fr

Tania Richmond
email: tania.richmond@univ-tln.fr

Improved Timing Attacks against the Secret Permutation in the McEliece PKC

Keywords: code-based cryptography, McEliece PKC, side-channel attacks, timing attack, extended Euclidean algorithm

In this paper, we detail two side-channel attacks against the McEliece public-key cryptosystem. They exploit timing differences on the Patterson algorithm used for binary Goppa codes decoding in order to reveal one part of the private key: the support permutation. The first one is improving two existing timing attacks and uses the correlation between two different steps of the decoding algorithm. This improvement can be deployed on all error-vectors with Hamming weight smaller than a quarter of the minimum distance of the code. The second attack targets the evaluation of the error locator polynomial and succeeds on several different decoding algorithms. We also give an appropriate countermeasure.

Introduction

One of the main threats in the modern cryptography is the arrival of quantum computers. It was shown that cryptosystems based on large number factorisation could be compromised [START_REF] Shor | Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer[END_REF]. In reaction to this thread, several solutions have been proposed, such as hash-based cryptography, code-based cryptography, lattice-based cryptography, and multivariate cryptography. The new facts concerning the post-quantum cryptography are well discussed in [START_REF] Daniel | Post-Quantum Cryptography[END_REF].

Code-based cryptosystems were introduced in 1978 by Robert J. McEliece [START_REF] Mceliece | A public-key cryptosystem based on algebraic coding theory[END_REF]. Many variants were attacked and partially or totally broken. Up to now, none of the proposed variants seemed as strong and secure as the original McEliece public-key cryptosystem (PKC) using Goppa codes. Structural attacks managed to reveal the secret key and totally break variants that used the generalized Reed-Solomon codes [START_REF] Sidelnikov | On the insecurity of cryptosystems based on generalized Reed-Solomon codes[END_REF] or QC-LDPC codes [START_REF] Otmani | Cryptanalysis of a McEliece cryptosystem based on quasicyclic LDPC codes[END_REF].

As the Goppa codes still resist to structural attacks, they present a real interest in our approach. So we focus our attention on the cryptanalysis of the McEliece PKC using Goppa codes. More exactly on the side-channel attacks using time differences between two executions of the same task. The interest of timing attacks is both practical and theoretical: we avoid unsecured implementations and discover new attacks succeeding in a polynomial time. The main purpose of these types of attacks is to reveal a part of the secret key and a breaking point of an algorithm. The authors of such exploits usually end up by giving the necessary countermeasures and the secure variant of the algorithms.

In the case of the McEliece PKC using Goppa codes, most of the timing attacks were discovered since 2008. Falko Strenzke's articles mention several weak points mostly situated in the decoding algorithm [START_REF] Strenzke | Side channels in the McEliece PKC[END_REF][START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF][START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF]. Some of these can be repaired by an intelligent and cautious way of the programming manner such as proposed in [START_REF] Avanzi | Side-channel attacks on the McEliece and Niederreiter public-key cryptosystems[END_REF][START_REF] Strenzke | Side channels in the McEliece PKC[END_REF][START_REF] Daniel | McBits: fast constant-time code-based cryptography[END_REF].

All of the mentioned attacks were realised on a McEliece PKC implementation using the Patterson algorithm (cf. Fig. 1) for decoding Goppa codes. The number of error corrections in the Patterson algorithm is bounded: up to t errors can be corrected, where t is the degree of the Goppa polynomial.

Our contribution is to reveal a new timing attack against the error-locator polynomial (ELP) evaluation and to improve two existing attacks. In our new version of the two combined existing attacks, we detail how the relation between the two attacks is crucial in order to avoid eventual errors. The attacks are executed on the extended Euclidean algorithm (EEA) and exploit the number of iterations. As authors mentioned in [START_REF] Strenzke | Side channels in the McEliece PKC[END_REF][START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF][START_REF] Avanzi | Side-channel attacks on the McEliece and Niederreiter public-key cryptosystems[END_REF], the initial attacks are limited and may not allow the total break of the permutation. This limit is situated in the number of equations detected by their attack. We will use a new relation between the number of iterations in the two steps in order to expand the system and to fully determine the secret permutation. We will also give a single countermeasure, which is efficient to all types of attacks exploiting the EEA in this particular manner.

The second contribution is in giving a new timing attack against the ELP evaluation. The importance of this new attack is that it operates on the polynomial evaluation, applied in several decoding algorithms as the Patterson algorithm, Berlekamp-Massey algorithm or any general decoder for alternant codes. We will show that this attacks succeeds on several variants of the polynomial evaluation.

Background

Goppa codes

We will focus exclusively on binary Goppa codes in this paper, but it is easy to generalize our results to q-ary codes.

• Goppa polynomial:

g(x) is a polynomial over F 2 m [x] with deg(g) = t. • Goppa support: L = {α 0 , α 1 , .., α n-1 } subset of F 2 m s.t. g(α i ) = 0. • The syndrome polynomial associated to c ∈ F n 2 : S c (x) = n i=1 ci x+αi

Definition:

The Goppa code is:

Γ (L, g) = {c ∈ F n 2 | S c (x) ≡ 0 mod g(x)}
The syndrome polynomial S c (x) satisfies the following property:

S c (x) = ω(x) σ(x) mod g(x) σ(x) is called the error locator polynomial (ELP): σ(x) = t i=1 (x + a i )
, where ∀i ∈ {1, . . . , t}, a i ∈ L.

Properties: A Goppa code Γ (L, g) is a linear code over F 2 . Its length is given by n = |L|, its dimension is k ≥ n -mt, where t = deg(g) and its minimun distance is d ≥ t + 1.

Irreducible binary Goppa codes are defined by an irreducible Goppa polynomial g and admit the maximum length n = 2 m . We use this type of codes in the rest of the paper. The following notations are in correlation with these codes.

Notations:

In this paper, we will use the following notations:

• For the permutation of the support elements:

Π(L) = L = (Π(0), Π(1), . . . , Π(α i ), . . . , Π(α n-2 )
). • Let P (x) be a monic polynomial of degree t over F 2 m with t roots denoted a i :

P (x) = x t + S t t-1 x t-1 + S t t-2 x t-2 + . . . + S t 2 x 2 + S t 1 x + S t 0 ,
where the coefficients S i ∈ F m 2 are the usual symmetric functions:

S t t-1 = t i=1 a i , S t t-2 = t i=1 t j=1 j =i a i a j , . . ., S t 1 = t j=1 t i=1 i =j a i and S t 0 = t i=1
a i . The relations between coefficients will be exploited in Section 3.

The McEliece Cryptosystem

The McEliece PKC [START_REF] Mceliece | A public-key cryptosystem based on algebraic coding theory[END_REF] is composed by the three following algorithms.

Key generation:

The first step is to generate the support (the set of n = 2 m elements) and the Goppa polynomial g of degree t. Then, the parity check matrix can be built and brought to a systematic form: [I n-k |R] in order to recover a generator matrix G of the Goppa code. We randomly choose a non-singular k × k matrix S and a n × n permutation matrix Π, and compute the public k × n generator matrix G = SGΠ.

The key generation procedure outputs sk = (Γ (L, g), S, Π) and pk = (n, t, G).

Message encryption:

• Inputs: message m ∈ F k 2 , public key pk = (n, t, R T ).

• Output: ciphertext z ∈ F n 2 .
1. Randomly choose an n-bit error-vector with weight wt(e) = t; 2. Encode z = mG ⊕ e; 3. Return z.

Message decryption:

• Inputs: ciphertext z ∈ F n 2 , secret key sk = (Γ (L, g), S, Π).

• Output: message m ∈ F k 2 .
1. Compute z = zΠ -1 ; 2. Find m = mS from z ⊕ e using Decode(z ) with the secret code; 3. Compute m = m S -1 ; 4. Return m. Decode(.) is an alternant decoder.

Existing side-channel attacks:

There are several papers on side-channel attacks against the McEliece PKC and a quick review must be done in order to clear up the reader's understanding. Most of the attacks target the Patterson decoding algorithm and exploit several weaknesses as follow:

There are mainly two types of attacks classified by their goal:

1. Attacks recovering the secret message m [1,14,16]; 2. Attacks recovering (fully or partially) the secret key sk [START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF][START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Fast and secure root-finding for code-based cryptosystems[END_REF][START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF]. The attacks on steps and are able to determine some relations on the support elements by counting the number of iterations in the EEA. We improve them in Section 2. The attack on step reveals error positions using timing differences in the ELP evaluation. The attacker is able to find the error-vector with a certain non negligible probability. The basic idea is that two different

Step Ref. Countermeasure z = zΠ -1 S z (x) = H z (x t-1 , . . . , x 2 , x, 1) T S z (x) -1 mod g(x)
via EEA [START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF] control flow τ (x) = x + S z (x) -1 b(x)τ (x) ≡ a(x) mod g(x) [START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF][START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF] in EEA make sure deg(a)

≤ t 2 ; deg(b) ≤ t-1 2 deg(ri) = deg(ri-1) -1 via EEA and deg(τ ) = t -1 σ(x) = a 2 (x) + xb 2 (
x) e = (σ(α0), σ(α1), . . . , σ(αn-1)) ⊕ (1, . . . , 1) [START_REF] Avanzi | Side-channel attacks on the McEliece and Niederreiter public-key cryptosystems[END_REF][START_REF] Strenzke | Fast and secure root-finding for code-based cryptosystems[END_REF][START_REF] Strenzke | Side channels in the McEliece PKC[END_REF] polynomials, with some different degrees, are not evaluated in the same time. So the timing difference gives some information on the error-vector. We improve this attack in Section 3. In the rest of this paper, we assume that an attacker chooses a weight 0 < r < t for the error-vector e and we use the following notations: deg(g) = t and wt(e) = r.

2 Timing attack against double using of the EEA Goal: The attacker's goal is to recover the secret permutation Π.

Identification of a leakage:

The leakage is identified at steps and of the Patterson algorithm. This type of attack was already published in [START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF]. The two steps using the EEA are considered as independent parts. In this section, we propose to show the relation existing between both steps and thus attack them. Indeed, the main problem of previous attacks is the limited number of cases in which they can be exploited. They just can be applied on wt(e) ∈ {2, 4} as shown in [START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF] or wt(e) ∈ {2, 4, 6} as presented in [START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF]. The problem comes from a simple fact: the number of iterations is given by two conditions. One of the condition is that all quotients in the EEA must be polynomials of degree equal to one. So when this condition is not fulfilled the number of iterations could not be any longer controlled by the attacker. We will use N and N as notations for the number of iterations in the 3rd step (respectively 5th step) of the Patterson algorithm.

Motivations of our attack:

We show that using the relation between both steps allows us to completely control the number of iterations. The other contribution is in finding the relation between both steps and in using it for building a larger set of equations. We will show that we are able to extend the limited equation number of the system up to wt(e) = deg(g)

2 . The main interest is that instead of finding only equations involving the permutation of 2, 4 or maybe 6 elements, we can extend the search as much as necessary in order to discover the secret permutation. In terms of complexity, instead of enumerating all possible permutations, i.e. n! permutations, we reduce the complexity to the following expression:

p i=3 n i , where p ≤ deg(g) 2 -1
So for small values of p, we have:

p i=3 n i ≤ n p × (p -2) ≤ en p p × (p -2)
(where e = 2.718281828 . . . is the basis of the natural logarithms).

In order to make clearer the difference between the original attack and ours, we propose to give a lower bound for the complexity of the original attack:

n e n ≤ n! Finally: p i=3 n i ≤ en p p × (p -2) << n e n ≤ n!
So the original attack is exponential in the length of the code as a timing attack would only have a complexity in the maximum of the error-vector's weight needed for the attack (often extremely small in comparison with the code's length).

Scenario:

The attacker proceeds in the three following steps:

1. He chooses a random message m and computes c = mG; 2. He randomly chooses an error-vector e of small weight wt(e) < t (where t is the correction capacity of the code) and computes z = mG ⊕ e; 3. He sends z to an oracle (O), which outputs the message m and the number of iterations in steps and of the Patterson algorithm from Fig. 1.

Main idea: For wt(e) = 2p with p ∈ N, the attacker finds equations having the following form:

wt(e) i=1 Π(α i ) = 0. He will be able to build this type of equations with 0 < wt(e) < deg(g) 2 .

Conditions: The general assumption is that the attacker knows the public key pk, the order of all elements in the support L (L is supposed to be public, for example in the lexicographic order) and has access to an oracle O. These assumptions are the same as in previously mentionned works. We improve the attack in the same context. The oracle O is also able to give some extra informations: the timing for the whole algorithm or just for one particular step. We assume that the attacker can violate the procedure by adding wt(e) < t errors. The attacker is able to choose the number and the positions of errors.

Step in the Patterson algorithm

It was shown in [START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF] that the syndrome inversion leaks some information. The attack is based on the number of iterations used in the EEA, in order to compute the inverse of the syndrome polynomial S(x) modulo the Goppa polynomial g(x). It uses the following properties:

N ≤ deg(σ) + deg(σ ), for wt(e) < deg(g) 2 .
We will not detail here all conditions, as they are well explained in [START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF], but we only give some important facts in order to make things clearer and to prepare the attack. Let us consider the ELP

σ(x) = x r + S r r-1 x r-1 + S r r-2 x r-2 + • • • + S r 2 x 2 + S r 1 x + S r 0 , with r ≡ 0 mod 2. Then σ (x) = S r r-1 x r-2 + • • • + S r 3 x 2 + S r 1 .
In this case, the maximum number of iterations is given by the coefficient S r r-1 = r i=1 a i . So if S r r-1 = 0, we obtain N = 2r -2 and all quotients have a degree equal to 1. If S r r-1 = 0 and S r r-3 = 0, then N = 2r -4 and all quotients have a degree equal to 1.

Step in the Patterson algorithm

Locate the leakage: Two observations have to be done in order to understand and to locate the leakage point. The first one is about the number of iterations. It was proven (in [START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF]) that this number is (with a high probability):

N = N i=1 deg(q i ) = deg(b).
In the following paragraph, we will give some relations between τ (x), b(x), a(x) and σ(x) (given in steps , and in Fig. 1). We will prove some new relations. The new relations between these polynomials allow us to build the attack in such a manner that previous ambiguous cases are eliminated. These relations are crucial for better understanding of the entire decryption algorithm as they influence each step of the process and each particular form of the involved polynomials.

Properties: There are few properties useful for our approach:

1. If r ≡ 0 mod 2, then deg(a) = r 2 see [11]. 2. If r ≡ 1 mod 2, then deg(b) = r-1 2 see [11]. 3. If deg(τ ) ≤ r 2 ,

Number of iterations

We can see that the number of iterations in the EEA equals to deg(b), so we will focus on the form of the polynomial b(x), more exactly when r is even. Let:

σ(x) = x 2p + S 2p 2p-1 x 2p-1 + S 2p 2p-2 x 2p-2 + S 2p 2p-3 x 2p-3 + . . . + S 2p 2 x 2 + S 2p 1 x + S 2p 0 .
We separate odd powers from even ones and get:

σ(x) = (x 2p + S 2p 2p-2 x 2p-2 + . . . + S 2p 2 x 2 + S 2p 0 ) +(S 2p 2p-1 x 2p-1 + S 2p 2p-3 x 2p-3 + . . . + S 2p 1 x) σ(x) = (x p + S 2p 2p-2 x p-1 + . . . + S 2p 2 x + S 2p 0 ) 2 +x( S 2p 2p-1 x p-1 + . . . + S 2p 1 b(x) ) 2
So deg(b) is given by the coefficients S 2p 2i-1 with i ∈ {1, 2, . . . , p}. Therefore the number of iterations could be given by the same coefficients under an extra condition: all quotients have a degree equals to one. So we can distinguish p -1 possible cases depending on the coefficients, if their degree is equal to 1 in each iteration. Therefore:

         N = p -1 if S 2p 2p-1 = 0 N = p -2 if S 2p 2p-1 = 0 and S 2p 2p-3 = 0 N = p -3 if S 2p
2p-1 = 0 S 2p 2p-3 = 0 and S 2p 2p-5 = 0 . . . In all cases, the same assumption is made: the degree of the quotient is equal to 1 in each iteration. It means that we might have the number of iterations without any condition on coefficients.

Attack against the pair (N , N )

How it works. In this paragraph, we will explain how our attack works. We start by presenting the general relation for the pair (N , N ). Using Subsections 2.1 and 2.2, we get the following property:

General property: Let wt(e) = 2p < t/2. N , N (N , N ) = (4p -4, p -2) ⇒ 2p i=1 a i = 0
with probability P success . In Appendix C, a toy example is presented for a better understanding.

Success probability

The success probability P success is described by the following event: {All quotients have a degree equals to one}. If we consider all elements of our support as uniformly distributed variables and the independence of each step inside the EEA, under the initial assumptions, we have:

P success = P({N = 4p -4} ∩ {N = p -2}) = P({N = 4p -4})P({N = p -2}) = (1 -1 n ) N +N .
Experimental results show that for n = 2048 and wt(e) = 4, in order to find equations, have the following form: Π(α 1 ) + Π(α 2 ) + Π(α 3 ) + Π(α 4 ) = 0 the probability is equal to 0.998. It means that less that 0.2% of the cases are not exploitable among all possible cases under the condition: N = 4 and N = 0. In other words, each time that this combination is revealed, the probability of having a good equation for our attack is equal to 0.998 for the given parameters.

Experimental work

In order to validate the relations that we present in the previous paragraph for (N , N ), we used a Pari/GP implementation of the McEliece cryptosystem (the code will be publicly available). We computed a keypair, then encoded and decoded a given message multiple times, by checking the value of the couple (N , N ), searching for the valid combinations described above. We have also used different values for m, the extension degree of the finite field. We ran the algorithm until we got the specific combination about hundred times. Then, we obtained an average value for the necessary iterations required to get the searched combination. The results are presented in the following table: It means that for m = 7, we need to send to the oracle in average 127 different ciphertexts, in order to get the wanted relation (Π(α 1 ) + Π(α 2 ) + Π(α 3 ) + Π(α 4 ) = 0). In the case of the previous equation, the density equals in average to:

1 127 × P success = 1 127 × 1 - 127 128 4 .
Knowing one relation allows us, by fixing one of the positions, to reduce the number of ciphertexts that has to be sent to the oracle, it means that wanted relations are revealed more often as we progress in the attack. It also gives the first intuition on the structure of the permutation (see Appendix C). Attack implementation: In order to practically test our attack, we used the same software implementation. In order to reveal timings close to real values, we repeated the attack more than 10 6 times. We present the obtained results in the following table :   wt(e) Timings for expected attack equation (in sec.) Timings for random type equation (in sec.) Observation: We didn't give the timings for the odd values as they are constant and independant from the linear combinations between the permutations of the error positions. From Figure 3, we observe that there's a slight difference between the attack on this type of combinations and on randomly distributed combinations.

4
As we mentioned before for the random combinations, those with the maximum number of iterations are more likely to appear (the case when all coefficients are different from zero). So in this case, we have the timing difference required for our attack to succeed. In Subsection 2.7, we explain how the patch will work not only on this type of attacks but even on other types as the bit-flipping attacks.

Countermeasures

We have seen that it is possible to attack a system by knowing how many times the EEA is repeated. The number of iterations can go from 0 to t -1 in the syndrome inversion and from 0 to t/2 in the ELP determination. In order to avoid a correlation-finding from the number of iterations, we propose to introduce extra iterations into the EEA. The number of extra iterations should be chosen between 0 and a value that we call extra. The extra value is either t/2 or t -1, for the syndrome inversion or the ELP determination , respectively. The variable i contains the number of iterations realized in the first part of the secured EEA. We chose to use integer values in the extra EEA steps, in order to avoid divisions by zero that may occur if we keep the previous terms. The point is to keep computing things that are as computationally expensive as the original EEA, so that an attacker can't make the difference between true steps and extra steps. We present the proposal of the secured EEA in Figure 4.

The new security parameters: We recall the fact that normal security parameters do not take in consideration of timing attacks. Usually security parameters are given under the assumption of possible theorical attacks such as ISD [START_REF] Becker | Decoding random binary linear codes in 2 (n/20) : How 1+1=0 improves information set decoding[END_REF]

1. d ← d break 2. [b-1, b0] ← [0, 1] 3. [r-1, r0] ← [g(x), f (x)] 4. i ← 0 5. While deg(ri) > d do i ← i + 1 ri-2(x) = ri-1(x)qi(x) + ri(x) bi(x) ← bi-2(x) + qi(x)bi-1(x) end while 6. a(x) ← ri(x) 7. b(x) ← bi(x) 8. extra = f (t, d break ) 9. While i < extra do i ← i + 1 ri-2(x) = 3qi(x) + 5 bi(x) ← 5 + 6qi(x)
end while

Fig. 4. The secured extended Euclidean algorithm

For the first parameters, a timing attack with p = 6 would reveal a complexity less than 2 61 elementary operations, that is way lower than the security level of the original proposal. So for timing attacks, larger parameters have to be taken in consideration in order to maintain the same level of security. For example, in order to reach a 100 bit security level against this type of timing attacks, one should propose n = 131072 without countermeasure. The usual solution is not to increase values of parameters but to propose secure variant of the algorithm, variant that is not vulnerable to the specified attack. Our proposal is slower that the original algorithm, it operates (t -1) × O(1) for the syndrome inversion and t 2 × O(1) for the key equation (where O( 1) is the usual complexity for a division). Meanwhile, it is secure against timing attacks described below. The proof is very simple and it's based on the fact that this particular type of timing attacks are based on the number of iterations in the EEA. Since our algorithm performs the same number of iterations, no matter which relations are hidden between the polynomial coefficients, it can't reveal any of such secret relations. Once the countermeasure applied, we ran the same attack and got the following timings for selected Hamming weights (the average timings are presented for more than 10 7 simulations in Fig. 5).

wt(e) Timings for attack type equations (in sec.) Timings for a random type combination (in sec.) Interpretation: We observe that the protected implementation is impossible to attack (using the same techniques). We stress that the proposed countermeasure is also efficient in the case when an attacker wants to use previous techniques, like in [START_REF] Strenzke | A Timing Attack against the Secret Permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF][START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF].

Timing attack against the ELP evaluation

Goal: The attacker's goal is to find the private permutation Π.

Identification of a leakage:

A leakage is identified at step of the Patterson algorithm: the ELP evaluation. We recall that the ELP is denoted σ in Subsection 1.1. The attack is based on the fact that the form of the polynomial differs depending on the word to decode. We will prove that the algorithm's complexity is strongly related to the coefficients of σ(x). We will then perform a timing attack on the ELP evaluation and control the values of the coefficients of σ(x).

Motivations of our attack: One of the main motivations of our attack is that it can operate on all existing implementations of a general alternant decoder. It operates on the ELP evaluation, step that has to be computed in any decoding algorithm solving the so-called key-equation. We will give two basic algorithms for the ELP evaluation with some improvements and show that even with the published improvements our attack succeeds. We will choose the polynomial evaluation from right to left (the straightforward algorithm) and from left to right (the Ruffini-Horner scheme). Let the σ polynomial be of degree t. The first algorithm computes the result within 3t -1 operations (t additions and 2t -1 multiplications), whereas the second one does it within 2t operations (t additions and t multiplications). It was proven by V. Pan in 1966 [START_REF] Victor | On Methods of Computing the Values of Polynomials[END_REF] that the Ruffini-Horner's scheme [START_REF] William | A new method of solving numerical equations of all orders,by continuous approximation[END_REF] is optimal in terms of complexity. We will see that our attack works in the case of the first algorithm. We will give an improvement for a faster computation of this algorithm but still vulnerable to our attack. For the second algorithm, the attack is still successful with an extra condition: the attacker has to be able to detect whether at each step the algorithm computes the same number of operations or less. If this condition is fulfilled the attack works as well as in the first case. In the next paragraph, we will give more details about improvements that we used in our implementation.

The main idea of the improvement is to use the fact that some support elements have particular properties (e.g. 0 and 1). Knowing the fact that one coefficient equals zero speeds up the algorithm for operations like multiplication or sum has a fixed value if zero is taken as one of the input elements, the same thing happens within the multiplication by one. So we will exploit these properties in order to improve our implementation. Each time a coefficient equals one or zero it will be stored in a special table used afterwards by multiplication or addition. The case where a coefficient equals zero is rare and its probability has been studied in [START_REF] Dragoi | Polynomial structures in code-based cryptography[END_REF]. Nevertheless, each time there's a coefficient equal to zero we will no longer multiply it by the corresponding element because of the zero product. So we will use the predefined tables to get rid of the useless operations. We will proceed exactly the same way when the multiplication of an element has to be done when a coefficient equals to one. So each time we have one coefficient equal to zero, using our predefined tables, we get rid of two operations (one addition and one multiplication).

Scenario:

The attack scenario is the same as in the previous attack except for the last step. Indeed, the attacker gets the running time for the ELP evaluation in this section (step in Figure 1).

Idea: For wt(e) = 2, the attacker will find the positions of Π(0) and Π(1) (permutations of zero and one). After enough iterations, he will fix those two positions and repeat this attack with wt(e) = 3, he will then find the secret permutation Π (using exhaustive search for the remaining positions).

Conditions:

The assumptions are the same as in the previous attack except that the attacker does not know the order of the elements in the support L.

Success probability

As we said, in this attack we will only consider polynomials with a degree lower than three. For the case r = 3, we will give the full table of probabilities. We will start with the following general problem:

Problem: Let P (x) be a monic polynomial of degree r with r distinct roots over F 2 m . What is the probability that all its coefficients are different from zero?

Conclusion

In this article, we focused our attention on the cryptanalysis of the McEliece PKC using the binary Goppa codes. We showed the existing weak points in the Patterson decoding algorithm and determined the relations between the number of iterations in two different steps of this algorithm and the secret permutation. Since those relations were the main connection idea between the two extended Euclidean algorithm calls, we set up a timing attack based on this fact. The advantage of this attack is that it increased the probability of success by avoiding ambiguous cases, undetectable in previous attacks. The other advantage is that it allows higher expansion of the number of equations determined by the attacker in order to find the secret permutation.

The second important contribution of our article is a new attack that can be performed on several different decoding algorithms. It reveals that even intelligent variants of some polynomial evaluation algorithms might leak information and need to be patched or replaced. The ideas of these attacks might be reused in any further implementations (using the algorithms mentioned before), secure variants must be used in order to avoid any leakage point.

A Proof of the general properties

1. As deg(b) ≤ r-1 2 ⇒ deg(bτ ) ≤ r -1.
So there is no division by g(x) when searching for a(x). There are particular cases where zero iteration is performed in the ELP determination. Interpretation: A special case is when the errors have an even Hamming weight. Before the computation of the EEA algorithm, we could exactly give the degree of both polynomials (a(x) and b(x)) if τ (x) has a degree equal to half of the error-vector's weight. For example, if we compute τ (x) and find out that it has a degree equal to 4, we know exactly that an error-vector of Hamming weight equal to 8 was involved in the process. We also know that a(x) equals τ (x) at a constant close (because deg(b) = 0).

C Toy example for the EEA attack

Consider F 2 4 [x] = F2[x] x 4 +x+1 . The generator matrix G of the Goppa code and the support L = {0, 1, α, α 2 , . . . , α 14 } are public. Let m ∈ F k 2 be the message and O the decoding oracle. We notice that if L is public, one can find G(x) such that L = F2[x] G(x) . The other way is equaly true: if G(x) is public then one can easily find L. Suppose that the secret permutation is:

Π(L) = L = {α,
+ i2 + i3 = 0 0 if i1 + i2 + i3 = 0 • deg(b) = 1 if i1 + i2 + i3 = 0 0 if i1 + i2 + i3 = 0

E Proof of Answer from page 11

We have: S 3 2 = a 1 + a 2 + a 3 , S 3 1 = a 1 a 2 + a 1 a 3 + a 2 a 3 and S 3 0 = a 1 a 2 a 3 . We notice that: (0, α k , α j ) / ∈ E = {(α 1 , α 2 , α 3 ) ∈ L 3 | 3 i=1 α i = 0}. So:

P(S 3 2 = 0 ∩ S 3 1 = 0 ∩ S 3 0 = 0) = A 2 n-1 A 3 n .
In order to establish the full table of probabilities we need to compute:

P 1 = P(S 3 2 = 0 ∩ S 3 1 = 0 ∩ S 3 0 = 0).
This leads to the system of equations:

   a 1 + a 2 + a 3 = 0 a 1 a 2 + a 1 a 3 + a 2 a 3 = 0 a 1 a 2 a 3 = 0 ⇔    a 1 + a 2 = a 3 a 1 a 2 = (a 1 + a 2 )a 3 a 1 a 2 a 3 = 0 ⇔    a 1 + a 2 = a 3 a 1 a 2 = (a 1 + a 2 ) 2 a 1 a 2 a 3 = 0 .
Without loss of generality, we can set a 1 = α i and a 2 = α i+c , then:

   a 1 + a 2 = a 3 a 1 a 2 = (a 1 + a 2 ) 2 a 1 a 2 a 3 = 0 ⇔    a 3 = α i + α i+c α 2i+c = α 2i (1 + α 2c ) a 1 = α i , a 2 = α i+c ⇔ α 2c + α c + 1 = 0.
If m and 2 are coprime, the polynomial x 2 + x + 1 has no root in F 2 m , then P 1 = 0. Finally we obtain: P(S 3 i = 0) = P(S 3 i = 0 ∩ S 3 j = 0), ∀i = j ∈ {1, 2, 3}.

Fig. 1 .

 1 Fig. 1. Patterson algorithm: existing timing attacks and countermeasures

  then deg(a) = deg(τ ) + deg(b). 4. If deg(τ ) ≤ r 2 and deg(τ ) = 0, then wt(e) ≡ 0 mod 2. Proof. cf Appendix A. When wt(e) is odd: Let deg(g) be equal to 2p + 1, with p ∈ N. For an error-vector with Hamming weight wt(e) = 2k + 1, where k ≤ p -1, we have the following relations: deg(b) = k, deg(a) ≤ k and deg(τ ) ≥ 2p -k. See examples in Appendix B.1. When wt(e) is even: Let deg(g) be equal to 2p. For an error-vector with Hamming weight wt(e) = 2k, where k ≤ p, we have the following relations: deg(a) = k, deg(b) ≤ k -1 and deg(b) = 0 ⇔ deg(τ ) = k. See examples in Appendix B.2.

Fig. 2 .

 2 Fig. 2. Number of necessary iterations to get the combination for error-vectors of Hamming weight 4, 6 and 8

Fig. 5 .

 5 Fig. 5. Timings (in sec.) for decryption in the case of n = 2 11 and t = 10

2 .B. 1

 21 One knows that wt(e) ≡ 0 mod 2 ⇔ deg(a) ≥ deg(b) + 1. Suppose that deg(b) ≥ deg(a) ⇒ deg(τ ) = 0 ⇒ deg(b) = deg(a) ⇒ wt(e) ≡ 1 mod 2, what is in contradiction with assumption. B Examples for deg(g) = 20 With wt(e) ≡ 1 mod 2

  α 2 , α 3 , . . . , α 14 , 0, 1} = { i |i ∈ (1 . . . 16)} -1 st step: • The attacker asks O to decode all the z = mG ⊕ e with wt(e) = 1. N and N reveals the position of Π(0): 15 . • This is mainly due to: σ(x) = x. (We have τ (x) = 0 and S -1 (x) = x.) -2 nd step: • The attacker asks O to decode all the z = mG ⊕ e with wt(e) = 4. (The positions ( i1 , i2 , i3 ) are the three non-zero positions of e and i4 = 15 .) The couple (N , N ) = (4, 0) reveals all ( i1 , i2 , i3 ) such that i1 + i2 + i3 = 0. Here ( i1 , i2 , i3 ) ∈ {( 1 , 4 , 16 ), ( 3 , 14 , 16 ), . . . }. • deg(σ) = 4 and deg(ω) = 2 if i1

Fig. 6 .

 6 Fig. 6. Probability for r = 3

  . For example, for the McEliece PKC, the usual parameters are:

100-bit security n = 2048, t = 50, 128-bit security n = 2960, t = 56, 256-bit security n = 6624, t = 115. Inputs: f (x), g(x), d break and t. Outputs: a(x) and b(x) s.t. a(x) ≡ b(x)f (x) mod g(x)

-n-2 2 equations. To illustrate how this attack works, a toy example is given in Appendix D.

This problem was treated in [START_REF] Dragoi | Polynomial structures in code-based cryptography[END_REF] and the results show that the probability can be bounded. For the classical parameters of the McEliece PKC, i.e. n = 2048 and t ≤ 50, the authors obtain: P ≥ 0.95.

The case r = 3:

Answer : Let P (x) be a monic polynomial of degree 3 with three distinct roots over F 2 m and m = 1 mod 2. The probability P r=3 that all its coefficients are different from zero satisfies:

Proof. C.f. Appendix E.

Finding the permutation of the zero and one elements in the support

1. Consider the error-vectors e i with wt(e i ) = 1.

In this case, the error locator polynomial has the following form: σ(x) = x + a i , with a i ∈ L = {0, 1, α, . . . , α n-2 }.

If a i = 0, there is one addition (+) in the σ(x) evaluation. 2. Consider the error-vectors e i with wt(e i ) = 2.

In this case, the error locator polynomial has the following form:

We distinguish two possible cases:

if a i a j = 0 The case (b) leads to a computation of the polynomial evaluation with one extra addition (+) and the timings reveal all the couples (α i , 0). We can assume now that the position of Π(0) is known. 3. We fix this position and we seek for the position of Π [START_REF] Avanzi | Side-channel attacks on the McEliece and Niederreiter public-key cryptosystems[END_REF]. Since the polynomial σ(x) = x 2 + S 2 1 x, the fastest evaluation is obtained for the couple (Π(0), Π(1)) as there is only one addition (+) and one square computation.

Attack scenario when r = 3

We will consider error-vectors with Hamming weight that equals 3. The corresponding σ(x) polynomial has always one of the eight following representations:

if S 3 2 = 0 and S 3 1 = 0 and S 3 0 = 0 6. σ(x) = x 3 + S 3 1 x if S 3 1 = 0 and S 3 2 = 0 and S 3 0 = 0 7. σ(x) = x 3 + S 3 0 if S 3 0 = 0 and S 3 2 = 0 and S 3 1 = 0 8. σ(x) = x 3 if S 3 0 = 0 and S 3 2 = 0 and S 3 1 = 0 11 From Fig. 6 presented in appendix E, we obtain the four following cases:

Several cases can be eliminated by considering the fact that we accomplished the first step of the current attack, that is why we know the position of Π(0). If we consider all error-vectors where a i = 0 ∀i ∈ {1, 2, . . . , n -1} (i.e. 0 is not a root of P (x)), we reduce the possibilities for σ(x). The new form of the system is the following:

In all cases, x 3 must be computed so we will not consider this part in the timing differences. In the structure that computes the polynomial evaluation the fastest is the last one. But this case is performed only when S 3 2 = 0.

Finding the positions of two elements such that Π(α j )Π(α k ) = 1

In order to increase the number of equations in our system, we exploit the fact that (F 2 m ) * is cyclic. Recall: we know the positions of Π(0), Π(1) and Π(α 1 ) + Π(α 2 ) + Π(α 3 ) = 0. Without loss of generality, we choose to fix "Π(0)" on the first position and choose two other positions such that the sum is different from 1. We are able to do that because we know the position of "Π(1)" and the couples (α 1 , α 2 ) such that 1+α 1 +α 2 = 0. We get two new positions b 1 and b 2 such that b 1 + b 2 = 1. The error locator polynomial is:

This form is the fastest to be computed as there is one less multiplication compare to the other case.

System resolution

Number of equations:

We will give the number of linear and quadratic equations obtained by the attacker. Finding the positions of Π(0) and Π(1) reduces the search set to (n -2) elements.

• The first set of linear equations:

Equation type (1):

The last equation is determined by all other ones because for the last couple only one possible solution remains available. For instance, if the attacker finds ( n-2 2 -1) different equations, then the last equation can be directly determined.

• The second set of linear equations:

As the first set, the last one can be determined by all others. This comes from the fact that for the three positions, we fixed the position of Π( 1) as the first one. So we have (n -2) possibilities on the second position. But there are two repetitions for each (Π(1), Π(α j ), Π(α k ))-vector.

• The third set of quadratic equations:

The total number of equations for Π(α i )+Π(α j )+Π(α k ) = 0 including the second set is equal to (n-1)(n-2), as the third position is fixed and the two others are free and different. Here, the number of repetitions equals six. So we obtain (n-2)(n-1) -3 rd step:

• The attacker asks O to decode all the z = mG ⊕ e with wt(e) = 4. (The positions ( i1 , i2 , i3 , i4 ) are the four non-zero positions of e.)

The couple (N , N ) = (4, 0) reveals all ( i1 , i2 , i3 , i4 ) such that i1 + i2 + i3 + i4 = 0.

Here ( i1 , i2 , i3 , i4 ) ∈ {( 1 , 2 , 10 , 16 ), ( 2 , 3 , 13 , 16 ), . . . }.

• deg(σ) = 4 and deg(ω

• The attacker asks O to decode all the z = mG ⊕ e with wt(e) = 6. (The positions ( i1 , i2 , i3 , i4 , i5 ) are the five non-zero positions of e and i6 = 15 .) The couple (N , N ) = (8, 1) reveals all

• The attacker asks O to decode all the z = mG ⊕e with wt(e) = 6. (The positions ( i1 , i2 , i3 , i4 , i5 , i6 ) are the six non-zero positions of e.)

The couple (N , N ) = (8, 1) reveals all ( i1 , i2 , i3 , . . . , i6 ) such that i1

-Last step: The attacker has to solve the following system of quadratic equations in order to find the secret permutation:

Solving the system will allow to fully determine the secret permutation 3 , α 4 , . . . , 0, 1}.

D Toy example for the ELP evaluation attack

x 3 +x+1 . The Goppa polynomial G and the support L = {0, 1, α, α 2 , α 3 , α 4 , α 5 , α 6 } are public. Let m ∈ F k 2 be the message and O the decoding oracle. We notice that if L is public one can find G(x) such that L = F2[x] G(x) . The other way is equaly true: if G(x) is public then one can easily discover L. Suppose that the secret permutation is: Solving the system will allow to fully determine the secret permutation Π(L) = {α, α 3 , 1, α 4 , α 5 , 0, α 2 , α 6 }.


Gilles Madi-Wamba

Yunbo Li
email: yunbo.li@irisa.fr

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

Nicolas Beldiceanu
email: nicolas.beldiceanu@imt-atlantique.fr

Jean-Marc Menaud
email: jean-marc.menaud@imt-atlantique.fr

Cloud workload prediction and generation models

Cloud computing allows for elasticity as users can dynamically benefit from new virtual resources when their workload increases. Such a feature requires highly reactive resource provisioning mechanisms. In this paper, we propose two new workload prediction models, based on constraint programming and neural networks, that can be used for dynamic resource provisioning in Cloud environments. We also present two workload trace generators that can help to extend an experimental dataset in order to test more widely resource optimization heuristics. Our models are validated using real traces from a small Cloud provider. Both approaches are shown to be complimentary as neural networks give better prediction results, while constraint programming is more suitable for trace generation.

I. INTRODUCTION AND RELATED WORK

In the past ten years, the use of cloud computing technologies has significantly increased. This increase comes with a major problem, namely energy consumption. In 2010, the electricity used only by data centers accounted for between 1.1% and 1.5% of the world energy consumption. For the US, it even has increased to a level between 1.7% and 2.2% [START_REF] Koomey | Growth in Data Center Electricity Use 2005 to 2010[END_REF]. Up to 30% of the electricity consumption of a typical data center is dedicated to cooling the system [START_REF] Pelley | Understanding and abstracting total data center power[END_REF]. In [START_REF] Beloglazov | Energy-aware resource allocation heuristics for efficient management of data centers for cloud computing[END_REF], Beloglazov et al. showed that the electricity consumption of a data center is linked to its resources usage. A lot of research focuses on optimizing the resources usage of a data center. One common way to optimize resource utilization is to use workload prediction [START_REF]Multivariate Time Series ELM for Cloud Data Centre Workload Prediction[END_REF].

To predict the workload of a data center, in [START_REF] Ismaeel | Using ELM techniques to predict data centre VM requests[END_REF] and [START_REF]Multivariate Time Series ELM for Cloud Data Centre Workload Prediction[END_REF], Ismaeel et al. use a special type of neural network called an extreme learning machine. An extreme learning machine is a single-hidden layer feed-forward neural network, where the weights connecting inputs to hidden nodes are randomly assigned and never updated [START_REF] Huang | Extreme learning machine: theory and applications[END_REF]. In their work, as in many work in the literature, they assimilate the workload of a data center to the number of virtual machines running in the data center. Therefore, forecasting the data center workload amounts predicting the number of virtual machines requests. Our approach assimilates the workload of a data center to its CPU usage that we model with times-series. In our context this approach is more suitable as the specification of the machines of a data center directly involves their CPU capacities and consequently their energy consumption [START_REF] Beloglazov | Energy-aware resource allocation heuristics for efficient management of data centers for cloud computing[END_REF]. Moreover, it can help to implement energy-efficient techniques such as CPU over-commitment [START_REF] Beldiceanu | Towards energy-proportional clouds partially powered by renewable energy[END_REF].

Another issue encountered in designing resource optimization heuristics for data center is the availability of real workload traces to validate the proposed algorithms. In most of the cases, for reasons such as confidentiality, the real data are not available in enough quantities. One of the few available traces in the community comes from Google and provides a one-month trace for about 12,500 machines [START_REF] Google | Cluster data 2011[END_REF]. Yet, this kind of hyper-scale Clouds are only minor contributors to the global energy consumption of data centers (under 4%), while small-and medium-sized data centers account for half of this global electricity consumption [START_REF] Whitney | Data Center Efficiency Assessment[END_REF]. In this paper, we focus on these small and medium energy-hungry Clouds that present workload patterns different from the ones experienced in largescale public clouds [START_REF] Whitney | Data Center Efficiency Assessment[END_REF].

We present two complementary models, one based on neural networks and the other one based on constraint programming to address these issues: workload prediction and lack of real traces. Therefore, the contribution of this paper is twofold.

• To predict the CPU workload of a data center, we present and compare two machine learning models, which are respectively based on constraint programming and neural networks. To the best of our knowledge, constraint programming has not been used to predict a data center workload. • To extend any real workload data set, we present a workload traces generator. The generator uses the learned constraint model to generate time series that are similar to real workload traces. The rest of the paper is organized as follows: Section II presents the real workload traces that we are using in this paper. Section III presents the machine learning models, and Section IV gives the workload generator model. Section V concludes this work.

II. REAL WORKLOAD DESCRIPTION

Finding real traces for data center's activity is very difficult. For this, we have partnered with a french SME company: EasyVirt, a company specializing in virtualized data center analysis. In part of its activities, EasyVirt deploys software probes in the infrastructure of their clients and archives data in a MySQL database. These probes collect system resources' consumption of the physical servers and virtual machines. For confidentiality reasons, the companies in which the data were collected are not cited. Similarly, we can not distribute collected raw data. But, the analysis can be distributed. For this paper, we selected a representative trace of a mid-sized data center embedding 50 physicals servers for about 1,000 VMs.

More precisely, from a read-only account on the VMwareV-Center, the solution recovers static and dynamic information. The static information collected is: the Datacen-ter/Cluster/Server/VM architecture, static information of physical servers (server model, CPU model, CPU frequency, cores number, RAM available, etc.), static virtual machine information (vCPU (virtual CPU) number, allocated memory, reserved memory, VMware Tools status, VMDK size, etc). Following the recovery of this information, the probes dynamically monitor physical servers and virtual machines in order to get their consumption and virtual machines lifecycle. Various resources are monitored: processor, memory, network and disk. This monitoring is realized every 30 seconds, without impacting VMwareVCenter performance. Data are stored in a classical relational database (MySQL).

The selected trace represents six months of activity, for a database of about 2GB. All of the results presented below are based on these traces.

III. PREDICTION MODELS

In order to design a scheduling heuristic that aims at optimizing resource utilization, we need to know in advance the eventual distribution of the resource usage in time. To do so, we need a workload prediction model based on historical traces. This section presents two different approaches to build such a model. The first approach makes use of constraint programming while the second one uses a neural network.

A. Workload classification

As our dataset is consequent (six months, 50 physical servers), it can present various behaviors. As a preliminary step, we employ a clustering technique, widely used in data analysis, to group the traces into subsets of similar traces. K-means [START_REF] Patel | Impact of outlier removal and normalization approach in modified k-means clustering algorithm[END_REF] is one of the simplest unsupervised learning algorithm for data clustering.

Given a series (x 1 , x 2 , x 3 , ..., x n ), the parameter k(k ≤ n) is used for specifying the number of clusters to be created. The K-means algorithm is following the three steps:

1) Determine k centroid coordinates corresponding k clusters 2) Calculate the distance from each object to the k centroids 3) Assign each object to a unique cluster and ensure the sum of distance within a cluster is minimized. The objective of K-means [START_REF] Mohamad | Standardization and its effects on k-means clustering algorithm[END_REF] can be expressed as:

arg min s k i=1 xj ∈Si ||x j -µ i || 2 (1)
where, µ i is the mean of object in S i . Once this preliminary step has been done, one can start learning the characteristics of the traces in each subset.

B. Constraint programming model.

In [START_REF] Beldiceanu | Towards energy-proportional clouds partially powered by renewable energy[END_REF], we very briefly sketched a constraint programming model to predict the workload of a data center. This Section refines this model that will now be used both for prediction and for generating new time series workload (in Section IV). Given a set of real workload traces modeled as time series, we proceed as follows:

• Clustering step: we first clusterize them into p clusters, as already presented in Section III-A. This is motivated by the fact that the workload of a data center may both depend on the day (e.g. week, weekend) as well as of the type of services it provides at specific time periods. • Offline step: we extract from each cluster cl i , some key properties that correspond to typical features of the time series (e.g. its highest peak, its number of peaks). Using these key features, we build a model m(cl i ) for each cluster cl i . This offline preprocessing is the learning step.

Finally, the prediction model is given by ∪ p-1 i=0 m(cl i ), the union of all the models m(cl i ) of each cluster cl i .

• Online step: at any time t, the prediction model should be able to foretell in real time the evolution of the workload at time t + , where is a time lapse to determine. This rest of this section is structured in the following way:

(1) First we recall some background on time series, (2) then we present the offline step (3), and finally we describe the online step.

1) Background on time series: The references [3], [2] describe a large set of time-series constraints, where a time series constraint is characterized by the following concepts:

• Signature of a time series: the signature of a time series is a sequence of comparison operators taking values in the set {<, =, >}. Each element of the signature is obtained by comparing two adjacent input values.

• Pattern: a pattern is a regular expression over the alphabet {<, =, >}. To find a pattern occurrence in a time series, its signature has to be computed first. A pattern occurrence is a maximal occurrence of a sequence of characters from the signature that matches the regular expression of the pattern. Table I gives examples of patterns that may occur in a time series. Figure 1 provides a visual example of the peak pattern.

pattern regular expression increasing < increasing sequence < (< | =) * < | < increasing terrace < = + < summit (< |(< (= | <) * <)) (> |(> (= | >) * >)) plateau < = * > proper plateau < = + > strictly increasing sequence < + peak < (= | <) * (> | =) * > inflexion < (< | =) * > | > (> | =) * < steady = steady sequence = + zigzag (<>) + (< | <>) | (><) + (> | ><)
TABLE I: Examples of patterns and their corresponding regular expressions.

• Feature: given a pattern occurrence, a feature is a quantifiable property of the pattern.

• Aggregator: given one or more occurrence of a pattern p and a feature f of p, an aggregator is a function (e.g. min, max) applied to the different feature values of each occurrence of pattern p.

peak

< (= | <) * (> | =) * >
Fig. 1: Example of a peak pattern.

• Footprint: given a time series ts of length n and a pattern p, the footprint fp p (ts) of pattern p is a sequence of n values that identifies all occurrences of pattern p in the time series ts.

2) Offline step (learning step): Using constraint programming techniques, we analyze input workload traces to extract relevant features. For each cluster cl i a model m(cl i ) is built after a 3-step analysis of cl i . For each pattern p of interest and for each input time series ts belonging to the cluster cl i we proceed as follows:

(1) The number of occurrence of pattern p in the time series ts is computed. (2) The footprint f p p (ts) of the pattern p is computed.

(3) Different values for the aggregation of features of pattern p are computed.

At the end, all the results for all patterns of interest are put all together. They are analyzed to extract different ranges of variation for each characteristic. Among these ranges, we select a subset of relevant ones. Those relevant ranges constitute the model m(cl i ) of the cluster cl i .

Example III.1. Consider a cluster cl containing the three following times series of length 15 each: We consider the pattern peak , the feature width of the peak , and the aggregations sum and max . From steps (1) and ( 2) we obtain the following:

(nb_peak,2,[0,0,0,0,1,0,2,2,2,2,2,2,0,0,0]). (nb_peak,5,[0,0,0,0,1,0,2,0,3,0,4,0,5,0,0]). (nb_peak,4,[0,1,1,1,1,0,2,0,3,3,0,0,4,0,0]).

Each line of the form (nb pat, val , fp pat (ts i )). gives the number val of occurrence of the pattern pat in the time series ts i as well as the footprint fp pat (ts i ) of the time series ts i for the pattern pat.

From step (3) we obtain the following set of facts of the form (aggr feat pat, val , fp pat (ts i )), where each fact gives the value val of the aggregation aggr of different values of the feature feat for each occurrence of the pattern pat in the time series ts i .

(max_width_peak,6,[0,0,0,0,1,0,2,2,2,2,2,2,0,0,0]). (max_width_peak,1,[0,0,0,0,1,0,2,0,3,0,4,0,5,0,0]). (max_width_peak,4,[0,1,1,1,1,0,2,0,3,3,0,0,4,0,0]). (sum_width_peak,7,[0,0,0,0,1,0,2,2,2,2,2,2,0,0,0]). (sum_width_peak,5,[0,0,0,0,1,0,2,0,3,0,4,0,5,0,0]). (sum_width_peak,8,[0,1,1,1,1,0,2,0,3,3,0,0,4,0,0]).

We next put these results together and we obtain the following facts where each fact is of the form range(r , pat, prop, min, max ), where prop is a feature of the pattern pat. min (resp. max ) is the smallest (resp. largest) value that the feature prop takes when evaluated wrt each time series of the cluster cl . Finally r is the range obtained by subtracting min from max . range(3, peak,nb_peak,2,5). range (5, peak,max_width_peak,1,6). range (3, peak,sum_width_peak,5,8).

The range gives information on the amplitude of variation of a feature. Thus we can restrict the features of interest for the model of the cluster by setting a maximum allowed range.

Example III.2. If we set max range = 4 then the model m(cl ) of the cluster cl is:

m(cl) = { range(3, peak,nb_peak,2,5), range(3, peak,sum_width_peak,5,8).}
A time series ts is thus compatible with the model m(cl ) of the cluster cl if and only if :

• 2 ≤ nb peak (ts) ≤ 5 • 5 ≤ sum width peak (ts) ≤ 8
3) Online step, prediction step: This section presents how we use the model built in the offline step to predicts the next values of a time series at run time.

We first introduce the notion of prefix time series:

Prefix time series Given an index t < n (where n is the length of a time series) the prefix time series induced by t (pref (t)) is the time series x 0 x 1 . . . x t which represents the workload from time 0 to time t. At time t < n, each value x i (i ∈ [0, t]) is already known.

The prediction is done in three steps. Given a prefix time series (pref (t)), we first determine to which clusters it may belong. The second step gives an interval for the potential values of the prefix time series (pref (t)) at time t + . Finally, we refine the interval to make the prediction more precise.

(1) First we check the compatibility of the time series ts with each cluster model m(cl ). The time series ts is compatible with cluster cl i if the value for each characteristic of ts (e.g. number of occurrence of each pattern, footprint of each pattern, aggregation values) falls in the corresponding range of the model m(cl i ) of cluster cl i . (2) For each compatible cluster cl, we compute the interval of potential values at time (t + ), and note this interval by I cl (t + ). We make the union over the different compatible clusters and denote it by I(t + ). (3) We reduce the size of the interval I(t + ). To do this, we consider the center time series of each compatible cluster cl i and compute the footprint of the patterns strictly increasing sequence and strictly decreasing sequence to identify locations where these patterns occur. We use this information to reduce the interval I(t + ).

4) Evaluation of the constraint programming based model : For the evaluation of this model, we started from a data-set of 500 real world workload traces times series. We partitioned the dataset into 5 clusters (using the method presented at Section III-A) of respective cardinalities 142, 79, 142, 106, and 31. The model was constructed by learning from 70% of the time series of each cluster, and the prediction online step were done with the 30% remaining time series, that we call the test time series.

Those traces are of length 1008 each, meaning that each hour is represented by 42 time steps. To evaluate the cluster compatibility part of the prediction and the quality of the predicted interval, we reduced the resolution of the traces to one value per hour, that is each hour is represented by one timestep.

Each series is therefor of length 24 and is associated to a cluster. For each prefix of length k (with k ∈ [START_REF] Beldiceanu | Using finite transducers for describing and synthesising structural time-series constraints[END_REF]23]) of each test time series we perform the following benchmarks.

(1) The first benchmark presented in Figure 2 evaluates the percentage of cases were there is at least one cluster that is compatible with the prefix time series.

Percentage

Fig. 2: Percentage of cases were there is at least one cluster that is compatible with the prefix time series.

(

) 2
The second benchmark presented in Figure 3 evaluates the percentage of cases were the actual (k + ) th value of a test time series belongs to our predicted interval. From the curves depicted in Figures 2 and3, we observe that the quality of the prediction depends on the cluster compatibility part of the prediction as both curves have a similar evolution.

We measure the quality of our prediction using a conventional tool in data analysis, called RMSE that we now recall.

Definition The Root Mean Squared Error is a standard measure used to evaluate the estimation of an unknown value. Given a sample of k predictions pred (1 ), pred (2 ) . . . , pred (k ), and the a set of k observed values val (1 ), val (2 ) . . . , val (k ), the RMSE of the prediction is estimated by :

RM SE = 1 k k i=1 pred (k ) -val (k ) 2
To compute the RMSE of our prediction using the constraint programming based model, we considered the center of the predicted interval I(t + 1) to be the predicted value. Figure 4 presents the RMSE of the prediction for each prefix length from 3 to 23.

RMSE value

Fig. 4: Evaluating the RMSE of the prediction.

The average value of the RMSE for all prefix lengths is about 2800. This is satisfactory since the values of the time series used for these benchmarks range from 3000 to 29000. However the quality of the prediction relies on the quality of the classification of time series into clusters. Meaning that, the results may go worse if new time series that we fit into the system are not correctly classified. Since we are building a real time system, we can not neglect this aspect as new time series are also classified in real time and there is no new evaluation of the classification.

To overcome these issues, we elaborated a new model based on neural networks, the model is presented in Section III-C and does not need a clustering phase.

C. Neural networks model

Before describing our neural networks model, we recall some background on artificial neural networks.

1) Background on artificial neural network: An artificial neural network is a structured and interconnected group of nodes that reads an input in and computes an output. Figure 5 shows a simple neural network.

The Neural network of Figure 5 is structured as follows :

• One input layer having five nodes.

• One single hidden layer having three nodes.

• One output layer having a single node.

A neural network may have more than one hidden layer, and the number of nodes of each layer may vary. Each arrow of the neural network carries a weight w, and each node is associated with a value b called bias. Further, a neural network is associated with a function called the activate function. Figure 6 shows parameters associated with a neural network.

Input layer

Hidden layer

Output layer

x 1

x 2

x 3

x 4

x 5

Output

Fig. 5: A simple neural network. Using its weights, biases and activate function, a neural network computes the output y of a given input x = (x 1 , x 2 . . . x k ) with the formula y = f ( k i=1 w i x i + b j ), where w i is the weight of the edge from input i to the node j of the hidden layer, and b j is the bias of node j of the hidden layer. This formula assumes that the neural network has one single hidden layer, but can be adapted to more than one layer.

In general, we denote by W the k by m matrix of weights, by B the m by 1 matrix of biases, and by X the input vector of size k. This simplifies the expression to : y = f (W X + B) (see [START_REF] Nielsen | Neural Networks and Deep Learning[END_REF]). The commonly used activate function is the sigmoid function σ = 1 e -z where z = W X + B. In all what follows, we set the activate function f of our network to be the sigmoid function σ. The learning problem of our network consists in using the training data to find values for the matrices W and B in a way that minimizes the cost function C. To solve this problem, we use the gradient descent algorithm. Next section presents the gradient descent algorithm.

The gradient descent algorithm: We recall that in a neural network, the learning problem is to find weights W and biases B that minimizes a cost function of W and B. In our case the cost function is the mean squared error function c(W, B).

To find W and B that minimizes the cost function, we use the gradient descent algorithm [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]. In all what follows, we consider c to be a function of w and b.

The gradient descent is a step by step algorithm. At the beginning, w and b are initialized to some values. The initial values may be chosen randomly. Lets denote w 0 and b 0 the values of w and b at step 0. The next step is to update values of w 0 and b 0 to w 1 and b 1 such that c(w 1 , b 1 ) ≤ c(w 0 , b 0 ). The change from c(w 0 , b 0 ) to c(w 1 , b 1 ) is ∆c and is given by the formula ∆c = ∆c δw ∆w + ∆c δb ∆c. Let ∆v = (∆w, ∆c) and ∇c = ( ∆c δw , ∆c δb ) be the gradient of c, we have ∆c = ∇c.∆v.

We need to choose ∆v such that ∆c is negative i.e such that c decreases.

Let ∆v = -η∇c where η is a positive number, then ∆c = -η ∇c 2 . Since ∇c 2 and η are positive, then ∆c is negative.

The value η determines how high (or how small) is ∆c, and is called the learn rate.

Next section presents how we use such a neural network to learn a model from a set of input time series.

2) Learning time series with a neural network : This section presents the design of our neural network as well as how it is trained to predict future values of time series.

The core idea behind the training of a neural network that predicts future values of time series is to learn a mapping from a prefix time series to its next value.

Given a set of real data called training inputs of the form (X = pref (t), a(X)) where pref (t) is a time series of length t and a(X) is the value of that time series at the next time step t + 1, the trained network should be able to predict a value y(X) that is close to a(X). Formally, this is done by solving a minimization problem. Let C(W, B) = 1 2n X y(X) -a(X)

2 be a cost function where :

• n is the number of training input.

• W is the weight matrix.

• B the bias matrix.

• y(X) = f (W X + B) is the output of the network with activate function f given the input X. The function C is the mean squared error, also known as the quadratic error function. From its definition, we see that C becomes small when y(X) tends to a(X). So the objective is to find values of the matrices W and B that minimizes C.

Since the input to our neural network is a prefix time series, we designed our network to be compatible with the maximum length a prefix can take in our application. As we are dealing with time series of length 24 that represents the daily workload of a data center, we thus design our network with 24 input neurons. When we want to feed the network with a prefix time series of length t < 24, we append -1 to the prefix to obtain a time series of length 24. Also since we want our network to predict a single value, we design it to have a single output neuron.

To compute the gradient ∇c we used the backpropagation algorithm [START_REF] Alsmadi | Back propagation algorithm: the best algorithm among the multi-layer perceptron algorithm[END_REF], a commonly used algorithm in machine learning.

The next section presents the evaluation of our neural network based prediction model.

3) Evaluation of the neural network based model: For the evaluation of this model, we used the same data as in the case of the constraint programming model (Section III-B4). To train the network, we extract prefixes of length 2 to 23 from each time series of the learning set, and feed the network with couple comprising those prefixes together with the next value of the time series. The number of training examples is therefor increased, since for each time series of length 24, from the learning set we extract 22 prefixes of length from 2 to 23. The neural network has 24 neurons in the input layer, one neuron in the output layer and one single hidden layer. The number of neurons in the hidden layer as well as the learning rate η were experimentally set to 65 and 0.2 respectively.

RMSE value

Fig. 7: Evaluating the RMSE of the prediction with the neural network model.

Figure 7 presents the RMSE for each prefix length from 3 to 23. With this model, the quality of the prediction is twice better than the results in the case of the constraint programming model were the best case RMSE value is above 2000. While the neural network model performs better than the constraint programming model, Figures 4 and7, shows that both prediction models follow the same patterns. The explanation is that the neural network learns the patterns from the learning time series. The peaks in the RMSE prediction curve correspond to time points where there are many possible candidate patterns. The valleys in the RMSE prediction curve correspond to time points where the current pattern of the time series has been identified. When the pattern is identified, the prediction is more accurate than when it is not yet identified. This observation remains true for the peaks and valleys observed in the RMSE prediction curve of the constraint programming model in Figure 4, that is why both curves follow the same patterns.

We recall that the purpose of this model is to predict at real time the workload of the data center. Real time means that the prediction should be fast enough to be exploited by the resource optimization heuristic. The time is therefore a very important parameter to take into consideration when dealing with such model. We also want to show that our model scales correctly. To evaluate the scalability of our model, we increased the resolution of the workload traces, from 24 values per hour to 41 values per hour. These benchmarks were performed on a computer running Mac OS 10.10.5 Yosemite, with 16GB of memory and an Intel core i7 processor at 2.93 GHz. From Figure 8, we can see that, with our model, the prediction time is less than 1ms even for time series of 41 values per hour, that is time series of length 984. This speed is due to the linear complexity complexity of forward propagation in a neural network. The number of multiplications needed to compute the output of the activation function of each neuron is linear with the number of neuron. The complexity is thus

O(k * 65) = O(k)
where k is the number of neuron of the input layer, in our case, k is the length of the time series. These result contrasts withe the time needed for the learning of that same neural network. These results led to Section IV on the use of the constraint programming model to generate new workload traces.

IV. WORKLOAD GENERATOR

Sections III-B and III-C presented and evaluated two prediction models. On the one hand, the results established that the model with neural network performs better than the constraint programming model, and is also well adapted for real time prediction as it requires no data classification.

On the other hand, the advantage of the constraint programming model is that it learns a set of constraints that characterizes the time series of each cluster. This very same model can then be used to generate time series that are compatible with a given cluster. Real data on server workload are in general not available in large quantities. Such a generator is very useful as it can generate as many data as needed. The generated data have the same patterns as the real workload data from the learning sets. The two models are thus complementary, as each one has its strengths. This section describes how new workload traces are generated in two steps from real workload data set. The first step builds a constraint satisfaction problem out of the available real workload traces, while the second step solves the problem to generate new traces. Before detailing the two steps, we recall some background on the notion of constraint satisfaction problem.

A. Background on constraint satisfaction problems Constraint satisfaction problem A constraint satisfaction problem (CSP) [START_REF] Rossi | Handbook of constraint programming[END_REF] is a triple P = (X, D, C) where :

• X is an tuple of n > 0 variables, X = x 0 , x 1 , . . . , x n . • D is a tuple of n > 0 domains, D = D 0 , D 1 , . . . , D n such that x i ∈ D i . • C is a tuple of k > 0 constraints C = C 0 , C 1 , ..., C k .
We now define the notion of constraint.

Constrain Given a CSP P = (X, D, C), a constraint C i ∈ C is a pair R i , S i where S i is a subset of X and R i is a relation over S i . The relation R i specifies the tuples of values forbidden among the Cartesian product of the domains of the variables in S i .

Solution of a CSP Given a CSP P = (X, D, C), a solution S of P is a tuple v 0 , v 1 , ..., v n such that v i ∈ D i and for every constraint C i = R i , S i ∈ C, the relation R i holds.

We now present how to construct a CSP to generate new workload traces.

B. Building the constraint satisfaction problem, and generating new traces

Lets P = (X, D, C) be the CSP that we need to construct in order to generate new workload traces. We define X, D and C as follows:

• [The set of variables X] We set the number n of variables of X to be the length of the real workload time series from the learning set. In our case n = 24.

i.e X = x 0 , x 1 , . . . , x 23 . To facilitate the reading, this model example considers only two characteristics of a single pattern. From the model we have the following information:

-All the time series of cluster cl have at least 2 and at most 5 peaks, thus the amplitude of the range of variation of the number of peaks is 3. -When summing the widths of all the peaks occurring in each time series of cluster cl, we have a value between 5 and 8. We translate this information into constraints, our problem thus comprises two constraints on the whole set of variables:

-

C 0 = R 0 , X with nb peak(R 0 , x 0 , x 1 , . . . , x 23 ) and R 0 ∈ [2, 5], -C 1 = R 1 , X
with sum width peak(R 1 , x 0 , x 1 , . . . , x 23 ) and R 1 ∈ [START_REF] Beloglazov | Energy-aware resource allocation heuristics for efficient management of data centers for cloud computing[END_REF][START_REF] Carlsson | SICStus Prolog user's manual[END_REF]. Using a constraint programming solver, we find the solutions S i = v i,0 , v i,0 , ..., v i,23 of the CSP C. Each solution S i is a time series of length 24 that has the same patterns learned from the real workload traces. The next section presents the evaluation of the quality of the time series generated by this model.

C. Evaluation of the CSP-based workload generator

To evaluate the quality of the time series generated, we ensure that it respects all the properties of the time series from the learning set. Those properties consist of all the features learned from the real workload traces. To solve the CSP model and generate solutions, we use a CSP Solver. For these benchmarks, we used Sicstus Prolog Solver [START_REF] Carlsson | SICStus Prolog user's manual[END_REF], on a computer running Mac OS 10.10.5 Yosemite, with 16G0 of memory and an Intel core i7 processor at 2.93 GHz.

By construction of our CSP model from the learned features, and by definition of a solution of a CSP, every solution i.e. every generated time series respects all features that the model learned from each cluster. We confirmed this theoretical expectation by checking weather each generated time series of length from 24 to 984 respects the properties extracted from the respective learning sets. For each length of time series, there are at least 40 learned properties, modeled as constraints.

Further, we evaluated the time needed by the CSP to generate new time series according to the length of the time series. Figure 10 shows that for any resolution from 1 to 41, the time needed to generate a time series of corresponding length, that is of length 24 up to 984 is less than one second. This good performance is explained by the paradigm of propagation constraint [START_REF] Bessiere | Constraint propagation[END_REF] that is used to solve the CSP model and generate solutions.

D. Neural network-based workload generator

We also used the neural network to generate time series. To do so, we start from a small prefix, i.e. a prefix of length As expected, none of the time series generated with the neural network respects all of the learned characteristics. Another disadvantage in trying to generate a new time series with the neural network is that the neural network needs to start with a non-empty small prefix that it completed to create a full time series. On the other hand, the CSP model can rapidly generate full times from scratch that verify all learned constraints.

V. CONCLUSION

In this paper, we have proposed two original workload prediction models for Cloud infrastructures. These two models, respectively based on constraint programming and neural networks, focus on predicting the CPU usage of physical servers in a Cloud data center. The predictions could then be exploited for designing energy-efficient resource allocation mechanisms like scheduling heuristics or over-commitment policies. We also provide an efficient trace generator based on constraint satisfaction problem and using a small amount of real traces. Such a generator can overcome availability issues of extensive real workload traces employed for optimization heuristics validation. While neural networks exhibit higher prediction capabilities, constraint programming techniques are more suitable for trace generation, thus making both techniques complementary. Our future work includes providing a website to access on-demand datasets produced by our generator using various real workloads that cannot be made directly publicly available.

Fig. 3 :

 3 Fig. 3: Percentage of cases were the actual (k + ) th value of a test time series belongs to the interval I(t + 1) predicted.

Fig. 6 :

 6 Fig. 6: Parameters of a neural network.

Fig. 8 :

 8 Fig. 8: Prediction time using the neural network.

Fig. 9 :

 9 Fig. 9: Time needed to train a neural network

  • [The set of domains D] Let m be the smallest value taken by the real workload time series from the learning set, at any time point; and let M be the largest value taken by the real workload time series from the learning set, at any time point. m and M are such that, for any time series ts of any cluster, and for any index t ≤ 23, the t th value ofts is included in [m, M ]. D = D 0 , D 1 , . . . , D 23 where D 0 = D 1 = • • • = D 23 = [m, M ].• [The set of constraints C] In Section III-B we presented the construction of a constraint programming model that captures all the features of the workload traces. We recall that a model m(cl) of a cluster cl is of the form:m(cl_i) = {range(3, peak,nb_peak,2,5), range(3, peak,sum_width_peak,5,8).}

Fig. 10 :

 10 Fig. 10: Time needed to generate times series from the CSP

  3, from every time series from the set of test time series; we predict the next value and repeat the process until we obtain a series of desired length. The process is described in Algorithm 1 Algorithm 1 Generate new time series with the neural networkL i ← Learning set T i ← T est set n i ← N eural network S i ← set of generated time series Desired length is N Each time series is of length N n i .train(L i )Train the neural network for st ∈ T i do p ← pref (st, 3) p is the prefix of length 3 of st while length(p) < N do y ← n.predict(p)

		Predict the next value
	p ← p y	concatenate y to p
	end while
	S i ← p
	end for

ACKNOWLEDGMENTS

This work has received a French state support granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference Nb. ANR-10-LABX-07-01.

The authors would like to thank the EasyVirt company for providing them real workload traces from Cloud providers.