Jean-Philippe Bernard

Emmanuel Frénod

Antoine Rousseau

Paralic confinement computations in coastal environment with interlocked areas

This paper is in the continuity of a work program, initiated in Frénod & Goubert [FG07], Frénod & Rousseau [FR13] and Bernard, Frénod & Rousseau [BFRss]. Its goal is to develop an approach of the paralic confinement usable from the modeling slant, before implementing it in numerical tools. More specifically, we here deal with the multiscale aspect of the confinement. If a paralic environment is separated into two (or more) connected areas, we will show that is possible to split the confinement problem into two related problems, one for each area. At the end of this paper, we will focus on the importance of the interface length between the two subdomains.

Introduction

Paralic confinement is one of the most pertinent parameters controlling the features of living species in paralic environments (i.e. environments such as lagoons, estuaries, bays, etc.). It was first introduced by Guélorget & Perthuisot [START_REF] Guélorget | Le domaine paralique. Expressions géologiques biologique, et économique du confinement[END_REF][START_REF] Guélorget | Le confinement, paramètre essentiel de la dynamique biologique du domaine paralique[END_REF] and it is linked with nutrient concentration of water in the paralic environments. It was widely discussed and tested in Guélorget [START_REF] Redois | Influence du confinement sur la répartition des foraminifères benthiques : exemples de l'estran d'une ria mésotidale de Bretagne méridionale[END_REF], Barnes [Bar94], Frénod & Goubert [START_REF] Frénod | A first step towards modelling confinement of paralic ecosystems[END_REF] and Tagliapietra et al. [START_REF] Tagliapietra | A review of terms and definitions to categorise estuaries, lagoons and associated environments[END_REF]. Its knowledge in a given paralic environment is an important factor for supporting decision of decision makers acting on the paralic environment. For instance, it can be used to help and choose the exact localization of shellfish farms, or to estimate the impact of the building of a new dyke or dam.

Since the recent works of Frénod & Goubert [START_REF] Frénod | A first step towards modelling confinement of paralic ecosystems[END_REF], Frénod & Rousseau [START_REF] Frénod | Paralic Confinement: Models and Simulations[END_REF] and Bernard, Frénod & Rousseau [BFRss], we know that it is possible to develop a methodology to simulate numerically the confinement in any paralic environment. We now enter a phase of our work program which long term objective is to provide an operational tool to compute the paralic confinement in any point of any paralic environment on earth, only from bathymetry and oceanographic data.

Many questions need to be reached before achieving such an objective and we start here by tackling a modest (but important) one, which is related to the capability of computing separately paralic confinement in two connected areas of a given paralic environment. We shall particularly focus on the interface boundary condition that is required for such a coupling.

Modeling issues

Interlocked areas

Coastal environment is made of interlocked areas. For instance, if we take a look at a marsh in a Mediterranean lagoon, we face with the following cascade of areas: Atlantic Ocean -Gibraltar Strait -Mediterranean sea -lagoon entrance -the lagoon -marsh entrance -the marsh. Beside this, coastal environments may present a wide range of scales. For instance in the previously evoked cascade, the marsh is several tens of meters large, while the lagoon size is about ten kilometers. Those two scales are small when compared with the characteristic size of the Mediterranean sea, which is itself small with respect to the Atlantic Ocean dimension. Sizes of transition inlets -Gibraltar Strait and the lagoon entrance -need also to be taken into account.

As it will be recalled in the sequel (see also [START_REF] Frénod | Paralic Confinement: Models and Simulations[END_REF]), the numerical computation of the paralic confinement in coastal environment first requires the computation of the water flow essentially going from the ocean to the coastal environment far end, induced by the combined effect of evaporation, tide and fresh water inputs from the rivers. Once the flow is known, a tracer following this flow and undergoing diffusion, is then computed. At the end of the process, this tracer provides the value of the paralic confinement. Because of the wide range of scales appearing in coastal environments (see above), a confinement simulation may rely on several mathematical models that one needs to couple. In the case where the two coupled models are identical, we face a domain decomposition problem. Even if the coupling between the two (ore more) subdomains should actually be two-way, we will focus on one-way exchanges (from the open deep sea to the lagoon, and finally to the marsh). In other words, when we decompose a computational domain in two parts, we will consider a main lagoon (the part that is directly connected to the open sea) and a secondary lagoon (see Figure 1 below). The confinement field in the main lagoon will have to be computed accurately in a truncated domain Ω main = Ω \ Ω seg , while the simulation in the secondary lagoon Ω seg is nothing but another "classical" (i.e. monodomain) confinement simulation, with the main lagoon playing the role of the open deep sea. In the next subsection we recall the model previously developed to compute the paralic confinement field in a lagoon. In subsection 3.3, we introduce the domain decomposition Ω = Ω main ∪Ω seg and the related issues. We pay a particular attention to the interface (and the related boundary conditions) between the main and secondary lagoons.

Domain and equations for paralic confinement computation in a lagoon

We consider (see Figure 2) a lagoon that is a cylinder with base a regular, connected and bounded domain Ω ⊂ R 2 with boundary ∂Ω. This boundary is shared into Γ in and Γ 0 with Γ in ∩ Γ 0 = ∅. Any point in Ω is denoted (x, y). The lagoon seabed is described by a piecewise continuous function b : Ω -→ R + , where b(x, y) represents the bathymetry level at the horizontal position (x, y) ∈ Ω. The water altitude h is such that h > sup Ω {b}, exluding outcrops. In summary, the geometrical model of the lagoon writes:

Lagoon = (x, y, z), (x, y) ∈ Ω, b(x, y) < z < h .
(1) Simulations will be conducted with δ/r 1 = 5%, 10%, 15% and 20%. The instantaneous confinement is -at a given time t and a given position (x, y) -the amount of time the water particle located at the position (x, y) at time t has spent inside the lagoon water mass. It is related to nutrient concentration at position (x, y) and time t. It is the result of two phenomena. In the first place, evaporation process generates a flow from the ocean to the lagoon far end. Essentially, ocean waters can be seen as having a very high nutrient concentration and when those waters travel towards the lagoon, they meet living species that take off nutrients. Then along the travel, their nutrient concentration decreases. Secondly, the nutrient concentration undergoes diffusion because of water molecular and eddy viscosities, but more pregnantly, because of the chemical precesses involved in the dissolution phenomena. This diffusivity is small, when compared with the average consequence of the flow. Nevertheless, it is important in places where the velocity of the flow is small, especially in the lagoon far end. Consequently, in order to compute the instantaneous confinement, we use a passive tracer g t advected by the water velocity field u u u and undergoing diffusion. As shown in Frénod & Rousseau [START_REF] Frénod | Paralic Confinement: Models and Simulations[END_REF], this model is compatible with any lagoon geometry (shape and bathymetry), with the only restriction that intertidal zones and seabed outcrops are not taken into account. The idea developped in [BFRss] is to solve the following advection-diffusion problem: for any time t > 0 and given a sufficiently large time T , the solution g t = g t (τ, x, y) of

Ω main Ω seg Ω Γ in Γ δ Γ 0 Ω main Ω seg Ω Γ in Γ δ Γ 0
                     ∀0 < τ < T, ∀(x, y) ∈ Ω, ∂g t ∂τ (τ, x, y) + u u u(t -T + τ, x, y) • ∇g t (τ, x, y) -ν∆g t (τ, x, y) = 0, ∀0 < τ < T, ∀(x, y) ∈ Γ in , g t (τ, x, y) = T -τ, ∀(x, y) ∈ Ω, g t (0, x, y) = T, (2)
is such that g t (T, x, y) is a good approximation of the value of the instantaneous confinement at time t ∈ R + and position (x, y) ∈ Ω. Here τ is a variable related to the time to spend into the lagoon and ν is the small (when compared with the average value of u u u) diffusivity coefficient that models the nutrient hability to spread out the water. In system (2), the water velocity field u u u(t, x, y) may be induced by several phenomena (such as evaporation, tide, river input, etc.), which are modeled by the generic function θ. We may compute this field by solving the following equation:

       -∇ • (h -b)u u u (t, x, y) = θ(t, x, y), ∀t > 0, ∀(x, y) ∈ Ω, ∇ × u u u = 0, ∀t > 0, ∀(x, y) ∈ Ω, u u u • n = F in (t, x, y), ∀t > 0, ∀(x, y) ∈ Γ in , u u u • n = 0, ∀t > 0, ∀(x, y) ∈ Γ 0 , (3)
where n stands for the unitary vector orthogonal to ∂Ω pointing outside Ω and where F in is a function defined on Γ in such that :

Γ in (h -b)F in (t, x, y) dl = Ω θ(t, x, y) dx dy, ∀t > 0. (4)
We remark that the velocity field u u u can be separated in several "elementary" velocity fields, each of those being solely induced by one single process Consequently, depending on those processes, the function θ can model one or several phenomena. The sytem of equations (3) is solved thanks to its corresponding velocity potential formulation (see [START_REF] Frénod | Paralic Confinement: Models and Simulations[END_REF]). Provided that ∇ × u = 0, we write u = ∇ψ and solve the following Laplace equation for ψ:

           -∇ • (h -b)∇ψ = θ on Ω, ( 5a
) ∂ψ ∂n = F in on Γ in , (5b)
∂ψ ∂n = 0 on Γ 0 . (5c)
In subsection 3.3 and section 4 below we will consider equations for the velocity in truncated regions of the lagoon (e.g. equations ( 10) and ( 6)). Naturally, these equations will always be solved thanks to the potential formulation, even if this is not explicitly specified.

Domain decomposition for a lagoon with a secondary lagoon

When in a lagoon a clearly separated entity -so called secondary lagoon -exists, we want to split problem (2), (3) and (4) into two problems -a first one set in the secondary lagoon and another one set in the remainder of the lagoon -being connected by conditions on the secondary lagoon entrance. Naturally, we want the concatenation of the results to approximate a solution of the system (2), (3) and (4) set in the whole lagoon with a good accuracy. The way to account for this situation, in what concerns the geometrical aspects, consists in sharing lagoon Ω into three parts (see the left picture in Figure 2): the secondary lagoon Ω seg , the main part of the lagoon Ω main and their common boundary Γ. They are such that Ω seg and Ω main are open subsets of Ω, Γ = ∂Ω seg ∩ ∂Ω main and Ω seg ∪ Ω main ∪ Γ = Ω. Moreover we assume that the secondary lagoon is not near the lagoon entrance, which is translated by ∂Ω seg ∩Γ in = ∅.

The key problem -and the most difficult one -is the obtention of the water flow and of the tracer within the main part of the lagoon without computing them in the secondary lagoon. We now focus on this issue and, in a first place, we build the system of equations to compute the velocity of the water flow in the main part of the lagoon Ω main . The two first equations of problem (3) are retained, because they describe the physics of water transport. The third equality of (3) that translates that water cannot escape from the secondary lagoon through the shore will also be kept. Thirdly, a condition on the interface -that we will discuss hereafter -will be written. Hence, in the main part of the lagoon Ω main , we write the following problem that allows us to obtain the velocity

           -∇ • (h -b)u u u main (t, x, y) = θ(t, x, y), ∀t ∈ R, ∀(x, y) ∈ Ω main , ∇ × u u u main = 0, ∀t ∈ R, ∀(x, y) ∈ Ω main , u u u main • n = F in (t, x, y), ∀t ∈ R, ∀(x, y) ∈ Γ in , u u u main • n = 0, ∀t ∈ R, ∀(x, y) ∈ Γ 0 , u u u main • (-n trans ) = -F (t, x, y), ∀t ∈ R, ∀(x, y) ∈ Γ, (6)
where vector n and function F in have the same definitions as in system (3) and equality (4), where n trans stands for the unitary vector, orthogonal to Γ, and pointing inside Ω main (or outside Ω seg ). The function F is to be determined. It is clear that if the solution u u u of (3) were known, we would choose F = u u u • n trans and then we would obtain a solution u u u main of (10) that would be such that u u u main = u u u |Ω main . Yet, we work under the assumption that the solution of (10) is not known (we indeed want to compute u u u main to have the value of u u u in Ω main ). Anyway, using the Laplace-Neumann compatibility condition (the same that brings us to write (4) for system (3) for u u u), we know an information on F which is

Γ (h -b)F (t, x, y) dl = Ω seg θ(t, x, y) dx dy, ∀t ∈ R, (7)
and that translates that the quantity of water entering Ω seg -and so leaving Ω main -through Γ compensates for what is consumed by the process modeled by θ over Ω seg . Knowing this information, we can consider that the missing information is the profile (or the shape) of F along the interface Γ. Approximations of this profile, that are classical and known as giving proper results, can be used. For instance, we can chose F as being constant along Γ or being a Poiseuille profile (see [START_REF] Jlm Poiseuille | Le mouvement des liquides dans les tubes de petits diamètres[END_REF]).

Having u u u main on hand, and then considering that we consequently know u u u on Ω main with a good accuracy, in order to compute the passive tracer given by (2) only on Ω main , we will consider the following problem:

                       ∀0 < τ < T, ∀(x, y) ∈ Ω, ∂g main t ∂τ (τ, x, y) + u u u main (t -T + τ, x, y) • ∇g main t (τ, x, y) -ν∆g main t (τ, x, y) = 0, ∀0 < τ < T, ∀(x, y) ∈ Γ in , g main t (τ, x, y) = T -τ, ∀(x, y) ∈ Ω, g main t (0, x, y) = T, (8)
which straightforwardly comes from (2) replacing u u u by u u u main . This system has to be coupled with boundary conditions on the interface Γ. As the diffusivity coefficient ν is small, following Halpern [START_REF] Halpern | Artificial boundary conditions for the linear advection diffusion equation[END_REF], choosing the following Neumann condition

∂g main t ∂(-n trans ) (τ, x, y) = 0, ∀0 < τ < T, ∀(x, y) ∈ Γ, (9)
will give a solution g main t which will correctly approach g t over Ω main .

Once brought a way to tackle the key problem, we can notice that we can implement a way to obtain the water flow and the tracer within the secondary lagoon. As a decision was made concerning the profile of F , and so concerning F on Γ, we can write the following system to obtain the velocity field u u u seg in the secondary lagoon Ω seg :

       -∇ • (h -b)u u u seg (t, x, y) = θ(t, x, y), ∀t ∈ R, ∀(x, y) ∈ Ω seg , ∇ × u u u seg = 0, ∀t ∈ R, ∀(x, y) ∈ Ω seg , u u u seg • n = 0, ∀t ∈ R, ∀(x, y) ∈ Γ 0 , u u u seg • n trans = F (t, x, y), ∀t ∈ R, ∀(x, y) ∈ Γ, (10)
where n has the same definition as in system (3) and where n trans and F are the one set to solve system (6). On the other hand, g main t on interface Γ can be computed as a result of system (8) and gives the value of the tracer on this interface. This function can be used as a Dirichlet boundary condition for the problem giving g seg t , which will be close to g t in Ω seg . This problem reads:

                       ∀0 < τ < T, ∀(x, y) ∈ Ω, ∂g seg t ∂τ (τ, x, y) + u u u seg (t -T + τ, x, y) • ∇g seg t (τ, x, y) -ν∆g seg t (τ, x, y) = 0, ∀0 < τ < T, ∀(x, y) ∈ Γ in , g seg t (τ, x, y) = T -τ, ∀(x, y) ∈ Ω, g seg t (0, x, y) = T, (11)
with

g seg t (τ, x, y) = g main (τ, x, y), ∀0 < τ < T, ∀(x, y) ∈ Γ. (12)
4 Numerical simulations

In this section, we present numerical simulations of the lagoon described in Figure 1 (see Section 3). These simulations were performed with the finite element method implemented in the FreeFem++ software [START_REF] Hecht | FreeFem++ manual[END_REF]. As in [BFRss] the velocity equation is solved thanks to a Laplace equation on the velocity potential (u u u = ∇ψ and we use P 2 elements for ψ), whereas the advection-diffusion equation on g t is solved thanks to P 1 elements. We consider four different configurations, in order to enhance the importance of the interface width |Γ|. The lagoon Ω is equally split in two parts Ω main and Ω seg , which is the most general (unfavorable) case. Indeed, in cases where |Ω seg | < |Ω main |, the truncation error obtained in the numerical simulations is lower.1

For each of these configurations, we perform a numerical simulation of confinement in the whole domain Ω, thanks to equations (2) and (3). The corresponding numerical solutions will be considered as reference solutions and denoted (u u u ref , g ref t ).

Simulations without interface information

We now consider the numerical simulation of confinement in the truncated domain Ω main (see Figure 1), in which we look for u u u main , g main t solutions of systems ( 6) and (8). For the sake of simplicity, we consider a flat bottom and set h -b ≡ 1, so that the boundary function F is such that:

F (t, x, y) = f pr (t, x, y) × Ω seg θ(t, x, y) dx dy ∀t ∈ R, ∀(x, y) ∈ Γ, (13)
where f pr denotes the (unknown) profile of the velocity along the interface Γ and is such that Γ f pr (t, x, y)dσ = 1. In the simulations below, we use a Poiseuille profile for f pr .

We enumerate in Table 1 the L ∞ norm of the relative error between the reference solution g ref t restricted to the subdomain Ω main and g main t (computed in Ω main from ( 8)-( 9)). This error is defined by

sup Ω main g main t -g ref t g ref t .
We observe from the results in Table 1 that the longer the interface, the larger the error. This is due to the lack of information we have at the interface (in particular on the velocity profile, see Section 3 and discussion above).

Simulations with interface information

We now reproduce the numerical simulations of Section 4.1, but instead of choosing a Poiseuille profile for the function f pr , we use the exact profile provided by the knowledge of u u u ref on the interface2 :

u u u main • n trans = u u u ref (t, x, y) • n trans ∀t ∈ R, ∀(x, y) ∈ Γ. ( 14
)
Remark 1 We obviously have that Γ u u u ref (t, x, y) • n trans dσ = Ω seg θ(t, x, y) dx dy, which means that the lack of information on Γ in Equation (13) concerns the velocity profile f pr rather than its average amplitude over the interface.

Then, thanks to the well-posedness of the corresponding velocity equations, uniqueness immediately insures that u u u main = u u u ref |Ω main , that is to say the knowledge of the velocity profile along the interface Γ insures the complete knowledge of the velocity in Ω main . As above, we enumerate the corresponding errors in Table 2, and as expected the errors also depend (in the same manner) on the interface width, but are notably lower than those of 1. The figure only illustrates the results indicated in the table: the plot corresponds to the infinite relative error (column 3) as the function of the interface width (column 2). Furthermore, we can illustrate the quality of the Neumann boundary conditions (9) used for the confinement equation with regard to the confinement diffusivity. Table 3 illustrates that when the diffusivity is low, Neumann boundary conditions efficiently approximate the exact transparent boundary conditions, which was already proved in [START_REF] Halpern | Artificial boundary conditions for the linear advection diffusion equation[END_REF].

Diffusivity ν L ∞ relative error 1.10 -1 0.010308 5.10 -2 0.00323477 1.10 -2 0.00267697 5.10 -3 0.00195359 Table 3: Relative error between confinements g main t and g ref t as a function of the diffusivity ν, in configuration 4 (worst case, large interface width δ = 20r 1 /100). The error increases with the diffusivity. The figure only illustrates the results indicated in the table: the plot corresponds to the infinite relative error (column 2) as the function of the diffusivity (column 1).

Conclusion

In this paper we are interested in the truncation of computational domains in confinement models. This is a very important issue both for classical domain decomposition problems and for the numerical simulation of confinement in limited areas of large lagoons, which we consider here. As soon as the domain truncation is done, the most important question to address is the search for artificial boundary conditions at the new boundary. It is known that their nature strongly depends on the PDE model that drives the considered process. Starting from the confinement model introduced in [FR13], we introduced some boundary conditions in order to limit the numerical error induced by the domain truncation. The chosen confinement condition is a classical homogeneous Neumann boundary condition, which is known to be accurate for small diffusivity values (see [START_REF] Halpern | Artificial boundary conditions for the linear advection diffusion equation[END_REF]). The truncation error is actually mainly due to the lack of knowledge of the velocity profile across the artificial boundary, particularly in the case where the interface and/or the secondary lagoon are large. It would be very interesting to evaluate how some partial informations on this velocity profile (provided by measurements) would improve the corresponding truncation error. We leave this to subsequent studies.

  , Frisoni & Perthuisot [GFP83], Guélorget et al. [GGLP90], Ibrahim et al. [IGF + 85], Debenay, Perthuisot & Colleuil [DPC93], Redois & Debenay

Figure 1 :

 1 Figure 1: Computational domain Ω. The lower part (below the interface Γ) is denoted Ω main . The upper part is denoted Ω seg . Left: large interface (δ = 20r 1 /100), where r 1 is the radius of Ω main . Right: tinier interface (δ = 10r 1 /100).Simulations will be conducted with δ/r 1 = 5%, 10%, 15% and 20%.

Figure 2 :

 2 Figure 2: Left: lagoon geometry, including a secondary lagoon. Right: A section of the lagoon geometry over a line going from the lagoon entrance to the lagoon far end.

t

  and g ref t with the same diffusivity ν = 0.01 and four different interface widths. The error increases with the interface width. The figure only illustrates the results indicated in the table: the plot corresponds to the infinite relative error (column 3) as the function of the interface width (column 2).

Table 1 :

 1 Relative error between confinements g main

	Configuration Interface width δ/r 1 L ∞ relative error
	1	20%	0.0279855
	2	15%	0.0212144
	3	10%	0.0133008
	4	5%	0.00627107

Table 1 .

 1 Configuration Interface width δ/r 1 L ∞ relative error

	1	20%	0.00267697
	2	15%	0.00204638
	3	10%	0.00134943
	4	5%	0.000690316

Table 2 :

 2 Relative error between confinements g main t and g ref t with the same diffusivity ν = 0.01 and four different interface widths. The error increases with the interface width, and is notably lower than in Table

These simulations have been performed, but for the sake of clarity we only present here the case where|Ω seg | = |Ω main |.

We could also imagine obtaining this information from measurements.

Acknowledgments

The authors are very grateful to A. Fiandrino for fruitful discussions related to this article.


Ronan Lashermes

Marie Paindavoine
email: marie.paindavoine@ens-lyon.fr

Nadia El Mrabet
email: elmrabet@ai.univ-paris8.fr

Jacques J A Fournier

Louis Goubin

Practical validation of several fault attacks against the Miller algorithm

Keywords: Pairing, Miller algorithm, fault attacks, blinding countermeasures, EM fault injection

Pairing based cryptography (PBC) is touted as an efficient approach to address usability and privacy issues in the cyberspace. Like most cryptographic algorithms, PBC must be robust not only against theoretical cryptanalysis but also against practical physical attacks such as fault injections. The computation of the Tate pairing can be divided into two parts, the Miller Algorithm and the Final Exponentiation. In this paper, we describe practical implementations of fault attacks against the Miller Algorithm validating common fault models used against pairings. In the light of the implemented fault attacks, we show that some blinding techniques proposed to protect the algorithm against Side-Channels Analyses cannot be used as countermeasures against the implemented fault attacks.

I. INTRODUCTION

Pairing based cryptography (PBC) is touted as an efficient approach to address usability and privacy issues in the cyberspace. As most cryptographic algorithms, PBC must be robust not only against theoretical cryptanalysis but also against practical physical attacks like fault injections. PBC uses bilinear mappings (or pairings) to construct cryptographic schemes for applications such as Identity-Based Encryption (IBE) [START_REF] Boneh | Identity-Based Encryption from the Weil pairing[END_REF], anonymous IBE, Marie Paindavoine's research for this paper was done during her internship with the CEA Tech.

Nadia El Mrabet wishes to acknowledge support from French project ANR INS 2012 SIMPATIC. one round tripartite Diffie-Hellman key exchange or searchable encryption [START_REF] Dutta | Pairing-Based Cryptographic Protocols : A Survey[END_REF]. A pairing calculation consists of two major steps namely the Miller Algorithm (MA) followed by the Final Exponentiation (FE). An exhaustive literature is available on the choice of curves and associated parameters for secure and efficient PBC implementation [START_REF] Freeman | A taxonomy of pairingfriendly elliptic curves[END_REF]. Recent advances in solving the Discrete Logarithm Problem (DLP) in small characteristics [START_REF] Joux | A new index calculus algorithm with complexity l(1/4+o(1)) in very small characteristic[END_REF] have reinforced the importance of focusing mainly on fields with a big prime characteristic. Issues linked to the resistance of PBC implementations against side channels and fault attacks have also been covered in the literature [START_REF] Mrabet | Differential Power Analysis against the Miller Algorithm[END_REF], [START_REF] Whelan | Side channel analysis of practical pairing implementations: Which path is more secure?[END_REF]. Most of the proposed fault attack schemes against PBC are theoretical and focus on the Miller algorithm [START_REF] Page | A Fault Attack on Pairing-Based Cryptography[END_REF], [START_REF] Whelan | The Importance of the Final Exponentiation in Pairings when considering Fault Attacks[END_REF], [START_REF] Mrabet | What about Vulnerability to a Fault Attack of the Miller's algorithm During an Identity Based Protocol?[END_REF] for which associated countermeasures have been proposed.

In this paper we explain how fault attacks have been implemented in practice against the Miller Algorithm independently from the final exponentiation. In the light of the implemented attacks, we analyse the countermeasures proposed in the literature to show that using blinding methods to add randomness to the intermediate calculations are not always efficient to thwart our attack.

The paper is organized as follows. We first provide a quick background on pairings and an overview of the current state-of-art of fault attacks on the Miller algorithm. We describe the principle behind and the practical settings for the fault attacks that we implemented. We then show that some blinding countermeasures are inefficient against our fault model and consequently look at how to properly secure the Miller algorithm.

II. BACKGROUND ON PAIRINGS

A. The Tate and Ate pairings

A pairing e is a bilinear and non degenerate map such that e : G 1 ×G 2 → G 3 , with G 1 , G 2 and G 3 three cyclic groups of the same prime order r. Let p be a prime number, E be an elliptic curve over F p and r a prime divisor of card(E(F p )). Efficient algorithms for pairings are realized with G 1 , G 2 subgroups of an elliptic curve E(F p k ) with point at infinity P ∞ and G 3 is the subgroup of the r th roots of unity in F p k , where k is the smallest integer such that r divides (p k -1). The integer k is called the embedding degree of E with respect to r.

The first pairings used in cryptography were the Weil and the Tate pairings [START_REF] Menezes | Reducing elliptic curve logarithms to logarithms in a finite field[END_REF], [START_REF] Frey | A remark concerning m-divisibility and the discrete logarithm problem in the divisor class group of groups[END_REF] and they have become a very active field of research with uses in innovative protocols like IBE [START_REF]Identity Based Cryptography[END_REF]. The latest and more efficient implementations are the Ate [START_REF] Hess | The Eta Pairing Revisited[END_REF], twisted pairings [START_REF] Hess | The Eta Pairing Revisited[END_REF], improved with the notion of optimal pairings [START_REF] Vercauteren | Optimal pairings[END_REF] and pairing lattices [START_REF] Hess | Pairing lattices[END_REF]. These pairings (that we call "Tatelike" pairings) are based on the model of the reduced Tate pairing defined by

e T : E(F p )[r] × E(F p k )/rE(F p k ) → µ r ⊂ F * p k (P, Q) → f r,P (Q) p k -1 r
, where f r,P (Q) represents the Miller function computed using the Miller algorithm and defined by the divisor r(P ) -([r]P ) -(r -1)(P ∞ ) evaluated at the point Q.

The Ate pairing has been developed in order to accelerate the computation time of a pairing, this "Tatelike" pairing is the one implemented and used in our experiments. Let π p be the Frobenius endomorphism on a point given by π p (X, Y ) = (X p , Y p ). Let t be the trace of Frobenius (#E(F p ) = p + 1 -t). The Ate pairing is defined by

P ∈ E(F p )[r] ∩ ker(π p -[1]), Q ∈ E(F p k )[r] ∩ ker(π p -[p]), K = t -1 and ate(Q, P ) = f K,Q (P ) p k -1 r
. Since the Miller algorithm is the central step for the most popular pairings, several optimizations have been proposed [START_REF] Scott | On the efficient implementation of pairing-based protocols[END_REF]. Algorithm 1 shows the Miller algorithm using the denominator elimination optimization. In PBC, if a secret is used, then it is one of the inputs of the pairing, either the point P or the point Q. From Algorithm 1, f 1 shall be called the Miller variable, the Doubling step shall consist of the operations at lines 4 and 5 while the Addition step is represented by lines 8 and 9.

Algorithm 1: Miller algorithm for the Ate pairing

Data:

K = (K n . . . K 0 ) 2 , P ∈ G 1 and Q ∈ G 2 ; Result: f K,Q (P ) ∈ G 3 ; 1 T ← Q ; 2 f 1 ← 1 ; 3 for i = n -1 to 0 do 4 T ← [2]T ; 5 f 1 ←-f 1 2 × h 1 (P ) ; 6 h 1 (x)
is the equation of the tangent at the point T ;

7 if Ki = 1 then 8 T ← T + Q ; 9 f1 ←-f1 × h2(P );

B. Barreto-Naehrig curves

Given the diversity of pairing implementations, it is not possible to study fault attacks on all of them. The 128-bit security level was chosen as usually recommended. The best current implementations with this security use Barreto-Naehrig (BN) curves [START_REF] Barreto | Pairing-Friendly Elliptic Curves of Prime Order[END_REF]. Furthermore, the fault attack we describe here is independent from the elliptic curve parameters and from the coordinates used. As an example, we used the parameters specified in [START_REF] Beuchat | High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves[END_REF] to construct our curve and our fields.

A BN curve is defined by the parametrized values t(x), r(x) and p(x) such that x fully defines the curve, with r(x) and p(x) prime integers. In our case x = 0x3FC0100000000000 and F p 12 is constructed through the following tower extension:

F p 2 = F p [u]/(u 2 -β), F p 6 = F p 2 [v]/(v 3 -u), F p 12 = F p 6 [w]/(w 2 -v), with β = -5 is a quadratic non-residue in F p , u is a cubic non-residue in F p 2 ,
and v is a quadratic nonresidue in F p 6 .

The base for F p 12 as a vector space over F p 2 is then (1, w, w 2 , w 3 , w 4 , w 5 ) and the base for F p 12 over F p is (1, w, w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 , w 9 , w 10 , w 11 ), and for R ∈ F p 12 , we denote :

R = R 0 + R 1 w + R 2 w 2 + R 3 w 3 + R 4 w 4 + R 5 w 5 + R 6 w 6 + R 7 w 7 + R 8 w 8 + R 9 w 9 + R 10 w 10 + R 11 w 11 ,
with w 2 = v, w 6 = u and w 12 = β.

The BN curve E that we use as an example is defined by the equation Y 2 = X 3 + 5 and the twisted curve E is given by Y 2 = X 3 -5/u = X 3 -u.

III. FAULT ATTACKS AGAINST PBC

An overview of fault attacks against PBC is given in [START_REF] Mrabet | Fault attacks on pairing-based cryptography[END_REF]. We briefly present here those attacks using the notations specified in Algorithm 1.

One of the first fault attacks on PBC was described by Page and Vercauteren [START_REF] Page | A Fault Attack on Pairing-Based Cryptography[END_REF] against the Duursma-Lee algorithm (efficient for supersingular elliptic curves over a finite field of characteristic 3 [START_REF] Duursma | Tate pairing implementation for hyperelliptic curves y 2 = x p -x + d[END_REF]). Page and Vercauteren use a random fault that occurs into the loop counter i of the Miller loop. The consequence is an observable difference in the number of iterations during the Miller algorithm. They recover the secret (point P or Q) by computing the quotient of two pairing results (one correct and one faulty). The quotient cancels out terms not influenced by the fault. Using the quotient and the equations of the Duursma-Lee algorithm, they solve a system thanks to the special form of the results of the Miller algorithm.

Whelan and Scott adopt a similar approach but with a different fault model [START_REF] Whelan | The Importance of the Final Exponentiation in Pairings when considering Fault Attacks[END_REF]. They propose to inject a fault in the Miller variable f 1 itself. They add a random value to the function line evaluated during the last iteration of the Miller algorithm. The attack consists in analysing the quotient between a valid and a faulty result.

The attack of Page and Vercauteren was extended by El Mrabet to the more general case of the Miller algorithm for any elliptic curve [START_REF] Mrabet | What about Vulnerability to a Fault Attack of the Miller's algorithm During an Identity Based Protocol?[END_REF]: the injected fault modifies the number of iterations of the Miller loop. Two faulty results corresponding to the Miller loop stopped at two consecutive iterations are needed. The author describes the recovery of the secret point using the quotient of those two intermediate results and the equations of the Miller algorithm.

More recently in [START_REF] Bae | Instruction fault attack on the miller algorithm in a pairing-based cryptosystem[END_REF], the authors proposed a fault model where the addition step is skipped with a fault attack. The instruction skip allows to recover h 2 which leads to an attack very similar to the previous ones mathematically.

To sum up, two approaches have been proposed to recover the value h 1 : either by corrupting the data flow (data modification) or by modifying the control flow (loop skip, test skip). Once we have h 1 , the Miller algorithm can be reversed to find the secret point.

Only one of the two input points is considered secret and the other one is public. This setting is notably representative of the level of knowledge which an attacker would have in an IBE scheme [START_REF] Boneh | Identity-Based Encryption from the Weil pairing[END_REF].

IV. IMPLEMENTED FAULT ATTACKS & RESULTS

Implementing a fault attack against the Miller algorithm consists of two steps: first finding h 1 which is in turn used to calculate the secret point. In this section, we illustrate how h 1 was found using faults before detailing how the secret point in the Miller algorithm can be recovered mathematically knowing h 1 (from [START_REF] Mrabet | What about Vulnerability to a Fault Attack of the Miller's algorithm During an Identity Based Protocol?[END_REF]).

A. Experimental set-up

The attacker wants to find h 1 (P ) at a known iteration. We use Electromagnetic (EM) pulses to generate faults in a microcontroller computing an Ate pairing. The platform used is similar to the one detailed in [START_REF] Moro | Electromagnetic fault injection: Towards a fault model on a 32-bit microcontroller[END_REF], [START_REF] Dehbaoui | Electromagnetic Transient Faults Injection on a Hardware and a Software Implementations of AES[END_REF].

The target of our experiment is a Cortex-M3 32 -bit microcontroller running at 56 MHz on which we run our own implementation of an Ate pairing following the algorithms in [START_REF] Beuchat | High-Speed Software Implementation of the Optimal Ate Pairing over Barreto-Naehrig Curves[END_REF] and coded in C language. The chip has not been designed as a secure chip but it bears some basic sensors for monitoring power and clock glitches which trigger hardware interrupts. These sensors are active during the experiments and are able to detect the EM pulses if they are too powerful. The board is underpowered at 2.8 V instead of 3.3 V in order to increase the sensitivity of the chip to the EM pulses.

The EM pulse injection platform used is composed of a pulse generator, a handmade coil antenna of diameter 1 mm (seen on Fig. 1), an XYZ motorized stage allowing to precisely (±1 µm) position the targeted chip with respect to the antenna, all this within a Faraday cage and with a computer controlling all those devices. At a precise moment determined by the attacker, the generator sends a voltage pulse, between -210 V and 210 V, to the coil antenna with rising and falling times of 2 ns and a duration of 10 ns. The coil is precisely placed just above the targeted position on the circuit. The EM pulse reaches the Power Ground Network of the chip and alters the computation performed by the device.

In order to find the right parameters of the fault injection set-up, we control and communicate with the chip through Keil µvision [START_REF] Arm | [END_REF] and its UVSOCK library. This allows us to access various internal registers giving us insights of what is going on after a fault injection even if an interrupt has been raised (i.e. the glitch has been detected). Experiments consist in finding the right parameters: position of the coil with respect to the chip, pulse amplitude, and the instant of the pulse emission with respect to a trigger raised by the targeted chip when reaching the operations of interest.

First we have to determine the locations where faults can be injected without interrupting the execution flow of the code. To do so, we vary the XY position of the EM probe over the chip's surface and observe how the chip behaves under an EM pulse stress. Three characteristic behaviours have been observed, depending on the spatial XY location, as illustrated in Fig. 2 With this approach we managed to find 2 set-ups for fault models corresponding to the modification of the control flow and that of the data flow (for recovering h 1 ):

• Controlled add model: a controlled value is added to h 1 through a fault injection. "Controlled" means that the value should be manageable: not totally random but with only a limited uncertainty allowing to individually test all possibilities. This fault model typically allowed us to implement the fault attack described in [START_REF] Whelan | Side channel analysis of practical pairing implementations: Which path is more secure?[END_REF]. • Loop skip model: the last iteration of the Miller algorithm is skipped. This fault model suits the needs for the fault attacks described in [START_REF] Page | A Fault Attack on Pairing-Based Cryptography[END_REF] and [START_REF] Mrabet | What about Vulnerability to a Fault Attack of the Miller's algorithm During an Identity Based Protocol?[END_REF].

Both these models have been theoretically proposed for fault attacks on pairings and allow to get h 1 . Since our implementation is an Ate pairing, the last iteration of the Miller algorithm has only a doubling step and no addition step.

B. First fault attack model: controlled add

In this attack, the secret point is Q and P is known to the attacker. In this experiment, we inject a fault during an addition in F p . The latter operation requires a multiword addition algorithm on the 32-bit chip.

We target one word of the addition and vary the instant of the injection of the pulse.

During a campaign (with carefully chosen parameters, seen on Fig. 3) of 1014 pulses emitted by the EM bench, 373 (37%) pulses induced an interrupt and 172 (17%) undetected faults were created on the least significant word. No fault was created on the other words, apart from the carry propagation due to the faulted least significant word.

The values obtained are:

• Correct value: 0x6EBCDA28 • 0x6EBCDA29 • 0xB75E6E10 ≈ 0x6EBCDA28/2 • 0xF 75F 792B • 0xF F F F F F F F
It is difficult to exactly explain what the exact effect of the electromagnetic pulses is since we do not have inside knowledge of the inner workings of the chip. We guess that in this case, the pulse affects the bus when data are retrieved from the RAM. The synchronisation may be altered which leads to 'shifts' within a word, or to the use of the bus pre-charge value of 0xF F F F F F F F .

It is possible to use a fault on the modular addition to recover h 1 (P ) if the attacker knows the point P . In this case, he can target the evaluation of one of the coordinates (R 0 , R 3 or R 4 ) of h 1 (P ) during the last iteration of the Miller algorithm. For example the The attacker can now recover h 1 (P ) by injecting a known fault e on the modular addition giving t * 0 = t 0 +e. This fault is propagated onto R 0 :

R * 0 = t * 0 •Y P = (t 0 +e)•Y P = R 0 +e•Y P = R 0 +∆ R0 .
(1) Since e and Y P are known to the attacker, ∆ R0 is known as well. At the last iteration of the Miller algorithm, we have:

f K,Q (P ) = f 2 1 × h 1 (P ) (2)
If the attacker is able to inject a known fault ∆ R0 in h 1 (P ), he recovers

f K,Q (P ) * = f 2 1 × (h 1 (P ) + ∆ R0 ). (3)
As he knows f K,Q (P ) * , f K,Q (P ) and ∆ R0 , he can find h 1 (P ):

h 1 (P ) = f K,Q (P ) × ∆ R0 f K,Q (P ) * -f K,Q (P ) (4)
If P is the secret and Q is known it is possible to obtain the same result with a fault on the last operation computing R 3 which is a modular subtraction.

C. Second fault attack model: loop skip

The goal of this experiment is to skip one, or several, iterations of the Miller algorithm. The operations of interest are the ones which update the loop counter and decide whether to quit the loop or not. NOP instructions have been added into the code when the trigger is raised in order to have an easier synchronisation between our code and the EM pulse (there is an incompressible delay of 80 ns or about 5 instructions). Additionally, the NOP instructions also clear the pipeline which lead to a better reproducibility of the faults. Our tests without the NOP instructions made the synchronisation work a lot harder and when a fault is injected, the probability of repeating it (with all parameters fixed) is lower. In order to have a clear presentation, the NOP instructions have been kept through the experiments in this paper. Yet, it is possible to remove them at the cost of harder experiments (and/or better equipment).

When injecting EM pulses during the execution of operations on the loop counter, several behaviours have been observed:

• Hardware crash: the chip is not responding and often needs to be reprogrammed. • Software crash: the chip is in an incoherent state (e.g. after a branch to a random location). • Interrupt: a problem has been detected by the chip which stops the normal course of the program. • Computation error: the program has performed a mistake but the chip is unaware of it. • Normal execution: the program has been normally executed. We fine-tuned our settings to perform loop skips in a "computation error" scenario: to do so, we first used a complete Miller algorithm in which we have included several counters which act as "snow-steps", the value of the counter easily telling us whether that instruction has been executed or not. We used these counters to identify the correct time of injection while maintaining the pulse amplitude at -210 V thus creating the highest stress. At the right time of injection, we could only observe interrupts at first. Finally by varying slightly the time of injection and the stress similarly to what can be seen on Fig. 3, we were able to inject a fault which induces the "skipping" of the branch instruction.

In the Ate pairing, the last iteration is a tangent evaluation only:

f K,Q (P ) = f 2 1 × h 1 (P ). (5)
Thus if we skip the last iteration, we obtain

f K,Q (P ) * = f 1 . (6)
Finally, h 1 is simply

h 1 (P ) = f K,Q (P ) (f K,Q (P ) * ) 2 . ( 7
)
By simply changing the targeted iteration, we are in fact able to leave the loop at whatever iteration we choose by skipping the corresponding branch instruction. We confirmed this by skipping the last two iterations simply by raising the trigger exactly one iteration earlier.

If it is possible to exit the loop at whatever iteration the attacker desires, the best choice would be to execute only the first iteration of the Miller algorithm. In this case

f K,Q (P ) * = h 1 (P ), (8)
it becomes unnecessary to compute a correct result of the Miller algorithm. This trick may be useful since accessing to the correct value of the Miller algorithm is in practice a feat, due to the final exponentiation.

D. Recovering the secret point from h 1

Now that we have seen how the value of h 1 (P ) was recovered in practice, we are going to illustrate how the secret point was derived from the latter value. We know that Q ∈ E(F p k ) and P ∈ E(F p ). In practice however we use the degree 6 twisted curve for the representation of Q: we simplify the notation by denoting Q as the point (x q : y q ) ∈ E (F p 2 ) instead of (x q • w 2 : y q • w 3 ) ∈ E(F p k ). Additionally, the point Q (and therefore

T ) is represented in jacobian coordinates (X Q : Y Q : Z Q ) which map to the affine representation (X Q /Z 2 Q : Y Q /Z 3 Q ).
The attacker can recover h 1 (P ) for the Ate pairing with

h 1 (P ) = (3X 3 T -2Y 2 T ) • w 6 + 2Y T Z 3 T y p • w 3 -3X 2 T Z 2 T x p • w 4 , (9) h 1 (P ) = R 0 + R 3 • w 3 + R 4 • w 4 . ( 10
)
with R 0 , R 3 , R 4 ∈ F p 2 (since w 6 = u ∈ F p 2 ) recovered through identification and T = [i]Q for some i known to the attacker. For the Ate pairing R 0 , R 3 , R 4 provide a system in F p 2 :

R 0 = (3X 3 T -2Y 2 T ) • u, R 3 = 2Y T Z 3 T y p , R 4 = -3X 2 T Z 2 T x p .
First, if P (for the Ate pairing) is the secret and Q is known, P can trivially be obtained with this system since T = [i]Q is known to the attacker and the system is linear in P coordinates. That is why from now on we focus on the case where the secret point is Q while P is known. The solution is barely more complex. The system now yields the univariate polynomial

R 2 0 β • Z 12 T + 4 R 0 u λ 2 2 -9λ 3 3 • Z 6 T + 4λ 4 2 = 0 (11)
with λ 2 = R3 2yp and λ 3 = -R4 3xp . This polynomial can be solved on F p 2 providing the value of Z T .

Once we know Z T , we use it into the initial system to obtain X T and Y T . The points which do not lie on the curve are eliminated. Finally, the possibilities for Q =

[i -1 ]T are computed.
If it is possible to control the loop at which the Miller algorithm returns its result, one would prefer the result after the first iteration. In this case no correct execution is needed to extract h 1 since f 0,Q (P ) = 1. This case was practically achieved just by moving the trigger to the first iteration.

E. The Final Exponentiation

The attacker must be able to access the faulty Miller results in our fault attack, which can be difficult in practice because of the presence of the final exponentiation for BN curves.

Several strategies have been proposed to deal with this final exponentiation. A scan chain attack has been proposed in [START_REF] Mrabet | What about Vulnerability to a Fault Attack of the Miller's algorithm During an Identity Based Protocol?[END_REF] to directly access the result of the Miller Algorithm. More recently [START_REF] Lashermes | Inverting the final exponentiation of tate pairings on ordinary elliptic curves using faults[END_REF], a scheme has been proposed whereby fault attacks can be used to reverse the final exponentiation. By showing that the Miller algorithm is vulnerable in practice, we open the way to a practical fault attack on the whole pairing with double faults. Yet how to properly override the Final Exponentiation in conjunction with a fault attack on the Miller Algorithm remains an open problem which has to be further studied, but this is not within the scope of this paper.

V. ANALYSIS OF THE COUNTERMEASURES

Once we have tested and demonstrated the practical feasibility of fault attacks on the MA, we hereby analyse some blinding countermeasures proposed in the literature and find out the loop skip fault model forbids to use the blinding countermeasures proposed against side-channel analyses as countermeasures against fault attacks.

A. Coordinates randomization

As proposed in [START_REF] Kim | Side channel attacks and countermeasures on pairing based cryptosystems over binary fields[END_REF] against side-channels analyses, instead of executing the Miller Loop with Q = (X Q :

Y Q : Z Q ) ∈ E(F p 2 ) (represented in jacobian coordi- nates), we execute it with Q = (λ 2 X Q : λ 3 Y Q : λZ Q ),
where λ ∈ F p 2 is a random blinding value. We note h

(λ) 1 , h (λ) 2 , f (λ)
K,P the line evaluations and the output of the Miller algorithm when performed with this blinding. When we include the λ in the tangent equation, we find that it is possible to factor it. For the doubling step the equation becomes h = λ 9i+12 h 2 , where i is an integer related to the number of iterations. The value of i can be found, but since this value has no influence on the result of the attack we remove it. Thus the result of the Miller loop is f

(λ) K,Q = λ a • f K,Q for some integer a.
In order to perform the fault attack, we need two different executions of the Miller algorithm. Thus we use two different blinding values, one for each execution. But actually, this adds only one unknown into the system:

h 1 (P ) (λ) = f (λ1)(τ +1) K,Q (P ) f (λ2)(τ ) K,Q (
P ) 2 = λ a 1 • f (τ +1) K,Q (P ) λ b 2 • f (τ ) K,Q (P ) 2 = λ a 1 λ b 2 • h 1 (P ). (12)
We name L = the new unknown and hence have:

h 1 (P ) (λ) = L • ((3X 3 T -2Y 2 T ) • w 6 + 2Y T Z 3 T y p • w 3 -3X 2 T Z 2 T x p • w 4 ), (13)
h 1 (P ) (λ) = R (λ) 0 + R (λ) 3 w 3 + R (λ) 4 w 4 . (14)
By identification of the decomposition in F p 2 , we obtain the system

R (λ) 0 = L • R 0 , R (λ) 3 = L • R 3 , R (λ) 4 = L • R 4 . ( 15
)
To this system, we add the equation derived from the fact that T lies on the curve E:

Y 2 T = X 3 T + 5Z 6 T . (16)
The system [START_REF] Hess | Pairing lattices[END_REF] and the equation ( 16) can be solved for the Ate pairing. To be solved, the resulting system requires the computation of the Gröbner basis which provides an equation of degree 12 in Z T only.

B. Miller's variable blinding

In this countermeasure again side-channels analyses, the line evaluation value is multiplied by a random element L of F p 2 for all iterations as suggested by Scott in [START_REF] Scott | Computing the tate pairing[END_REF]. Thus, we have:

h (λ) 1 = L • h 1 . (17)
As can be seen immediately, this countermeasure can be bypassed in the exact same way as the previous one with the system (15) and equation [START_REF] Scott | On the efficient implementation of pairing-based protocols[END_REF].

C. Protecting the Miller algorithm

In order to protect efficiently the Miller algorithm, countermeasures should take into account the fault models. We have shown that the loop skip model allows to bypass the coordinates randomization and Miller's variable blinding countermeasures even if the attacker is not able to replay a mask already used. Two other blinding techniques have been proposed in the literature, namely the additive blinding and the multiplicative blinding which are effective if the mask values are properly changed at each execution.

In the additive blinding technique [START_REF] Page | A Fault Attack on Pairing-Based Cryptography[END_REF], a mask is added on one of the input points and a second pairing is computed to remove the mask: e(Q, P ) = e(Q, P + M )e(Q, M ) -1 if Q is the secret point.

Another technique, called the multiplicative blinding [START_REF] Page | A Fault Attack on Pairing-Based Cryptography[END_REF], moves the security concern away from the pairing computation. Let a be a random value modulo r and b such that a • b = 1 mod r. Then e(Q, P ) = e([a]Q, P ) b . It appears that even if the attacker is able to recover [a]Q, he cannot find the secret point Q thanks to the ECDLP. Now the security of the system relies on the secrecy of a and b.

The protection of the Miller algorithm against the loop skip fault attack can also be realized with the protection of the loop itself. It can be enough to securely ensure that the correct number of iterations have been computed, i.e. by verifying the counter or with timing monitoring for examples. Moreover adding 'classical' countermeasures causing de-synchronisations (random executions within the Miller Algorithm) would complicate the experimental set-up of the attacks proposed in this paper.

VI. CONCLUSION

In this paper we describe how two (theoretical) fault attack models proposed in the literature against the Miller algorithm of a pairing calculation have been validated in practice against an Ate pairing implemented on a 32-bit general purpose processor. To our best knowledge, this is the first paper that validates in practice the fault attacks described against the MA. Such attacks do not put at risk the entire pairing calculation so long as the final exponentiation is considered "fault-proof" but it nevertheless proves that most of the theoretical attack schemes described in the literature against the MA can be implemented. In the light of the results obtained, we review the blinding countermeasures proposed in the literature and show that some of them are inefficient. This work shows that those blinding countermeasures should not be used on their own to protect the Miller Loop against fault attacks. Future work shall consist in testing other attack routes described against the Miller algorithm and other countermeasures. A complete fault attack on a pairing taking into account the final exponentiation would be of interest both theoretically and experimentally.

10 h2

 10 (x) is the equation of the line (QT );

Fig. 1 .

 1 Fig. 1. EM bench: the chip and the probe.

  : either no faults have been detected (the white strips in the figure); or interrupts are generated (the red/dark spots on the figure); or undetected faults have been inserted (the blue/light spots on the figure).

Fig. 2 .

 2 Fig. 2. XY map: red (dark) corresponds to interrupts and blue (light) to undetected faults (the white stripes and dots are visual artefacts where no faults have been injected).

Fig. 3 .

 3 Fig. 3. Delay (horizontal) and stress (vertical) variations: red (dark) corresponds to interrupts and blue (light) to undetected faults.

Fig. 4 .

 4 Fig. 4. Assembly code for the loop test SUBS r4,r4,#1 //index decrement CMP r4,#0x00 //loop test: is r4 == 0 BGE 0x080008C2 //>= branch to start

1 = λ 24i h 1 ,

 11 and for the addition step it becomes h (λ) 2