Olivier Bokanowski

Maurizio Falcone
email: falcone@mat.uniroma1.it

Smita Sahu
email: sahu@mat.uniroma1.

AN EFFICIENT FILTERED SCHEME FOR SOME FIRST ORDER HAMILTON-JACOBI-BELLMAN EQUATIONS *

Keywords: Hamilton-Jacobi equation, high-order schemes, -monotone scheme, viscosity solutions, error estimates

We introduce a new class of "filtered" schemes for some first order non-linear Hamilton-Jacobi-Bellman equations. The work follows recent ideas of Froese and Oberman (SIAM J. Numer. Anal., Vol 51, pp.423-444, 2013). The proposed schemes are not monotone but still satisfy some -monotone property. Convergence results and precise error estimates are given, of the order of √ ∆x where ∆x is the mesh size. The framework allows to construct finite difference discretizations that are easy to implement, high-order in the domains where the solution is smooth, and provably convergent, together with error estimates. Numerical tests on several examples are given to validate the approach, also showing how the filtered technique can be applied to stabilize an otherwise unstable high-order scheme.

1. Introduction. In this work, our aim is to develop high-order and convergent schemes for first order Hamilton-Jacobi (HJ) equations of the following form

∂ t v + H(x, ∇v) = 0, (t, x) ∈ [0, T ] × R d (1.1) v(0, x) = v 0 (x), x ∈ R d . (1.2)
Basic assumptions on the Hamiltonian H and the initial data v 0 will be introduced in the next section. It is well known that, in the one dimensional case, there is a strong link between Hamilton-Jacobi equations and scalar conservation laws. Namely, the viscosity solution of the evolutive HJ equation is the primitive of the entropy solution of the corresponding hyperbolic conservation law with the same hamiltonian. There are several schemes developed for hyperbolic conservation law (see references [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] [16], [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF], [START_REF] Gottlieb | Total variation diminishing runge-kutta schemes[END_REF]). Most of the numerical ideas to solve hyperbolic conservation law can be extended to HJ equations. Well known high-order essentially non-oscillatory (ENO) scheme have been introduced by A. Harten et al. in [START_REF] Harten | Uniformly high order essentially non-oscillatory schemes[END_REF] for hyperbolic conservation laws, and then extended to HJ equation by Osher and Shu [START_REF] Osher | High order essentially nonoscillatory schemes for hamilton-jacobi equations[END_REF]. ENO schemes have shown to have high-order accuracy although a precise convergence result is still missing and, for this property, they have been quite successful in many applications. Despite the fact that there is no convergence proof of ENO schemes towards the viscosity solution of (1.1) in the general case, convergence results may hold for related schemes, e.g. MUSCL schemes, as it has been proved by Lions and Souganidis in [START_REF] Lions | Convergence of MUSCL and filtered schemes for scalar conservation laws and Hamilton-Jacobi equations[END_REF]. Convergence results of some non monotone scheme have also been shown in particular cases [START_REF] Bokanowski | Convergence of a non-monotone scheme for Hamilton-Jacobi-Bellman equations with discontinuous initial data[END_REF]. Another interesting result has been proved by Fjordholm et al. [START_REF] Fjordholm | Arbitrarily high order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws[END_REF], they have shown that ENO interpolation is stable but the stability result is not sufficient to conclude total variation boundedness (TVB) of the ENO reconstruction procedure. In [START_REF] Fjordholm | ˙High-order accurate entropy stable numerical schemes for hyperbolic conservation laws[END_REF], a conjecture related to weak total variation property for ENO schemes is given. Finally, let us also mention [START_REF] Carlini | A weighted essentially nonoscillatory, large time-step scheme for Hamilton-Jacobi equations[END_REF], where weighted essentially non-oscillatory (WENO) schemes have been applied to HJ equations; the convergence proof of the scheme relies also on the work of Ferretti [START_REF] Ferretti | Convergence of semi-Lagrangian approximations to convex Hamilton-Jacobi equations under (very) large Courant numbers[END_REF] where higher than first order schemes are proposed in a semi-Lagriangian setting, yet with restrictive conditions on the mesh steps.

In this paper we give a very simple way to construct high-order schemes in a convergent framework. It is known (by Godunov's theorem) that a monotone scheme can be at most first order. Therefore it is needed to look for non-monotone schemes. The difficulty is then to combine non-monotonicity of the scheme and convergence towards the viscosity solution of (1.1), and also to obtain error estimates. In our approach we will adapt a general idea of Froese and Oberman [START_REF] Froese | Convergent filtered schemes for the Monge-Ampère partial differential equation[END_REF], that was presented for stationnary second order Hamilton-Jacobi equations and based on the use of a "filter" function. Here we focus mainly on the case of evolutive first order Hamilton-Jacobi equation (1.1), and an adaptation to the steady case will be also considered. We will design a slightly different filter function for which the filtered scheme is still an " -monotone" scheme (see Eq.2.12), but that improves the numerical results. Let us mention also the work [START_REF] Bokanowski | Value iteration convergence of -monotone schemes for stationary hamilton-jacobi equations[END_REF] for steady equations where some -monotone semi-Lagrangian schemes are studied.

The paper is organized as follows. In Section 2, we present the schemes and give main convergence results. Section 3 is devoted to the numerical tests on several academic examples to illustrate our approach in one and two-dimensional cases. A test on nonlinear steady equations , as well an evolutive "obstacle" HJ equation in the form of min(u t + H(x, u x ), u -g(x)) = 0 for a given function g are also included in this section. Finally, Section 4 contains our concluding remarks.

2. Definitions and main results. (2.1)

(A2) H : R d × R d → R d satisfies,
for some constant C ≥ 0, for all p, q, x, y ∈ R d :

|H(y, p) -H(x, p)| ≤ C(1 + |p|)|y -x|, (2.2)
and

|H(x, q) -H(x, p)| ≤ C(1 + |x|)|q -p|. (2.3)
Under assumptions (A1) and (A2) there exists a unique viscosity solution for (1.1) (see Ishii [START_REF] Ishii | Uniqueness of unbounded viscosity solution of Hamilton-Jacobi equations[END_REF]). Furthermore v is locally Lipschitz continuous on [0, T ] × R d .

For clarity of presentation we focus on the one-dimensional case and consider the following simplified problem:

v t + H(x, v x ) = 0, x ∈ R, t ∈ [0, T ],
(2.4)

v(0, x) = v 0 (x), x ∈ R. (2.5)
2.2. Construction of the filtered scheme. Let τ = ∆t > 0 be a time step (in the form of τ = T N for some N ≥ 1), and ∆x > 0 be a space step. A uniform mesh in time is defined by t n := nτ , n ∈ [0, . . . , N ], and in space by the nodes x j := j∆x, j ∈ Z.

The construction of a filtered scheme needs three ingredients:

• a monotone scheme, denoted S M • a high-order scheme, denoted S A • a bounded "filter" function, F : R → R.

The high-order scheme need not be convergent nor stable; the letter A stands for "arbitrary order", following [START_REF] Froese | Convergent filtered schemes for the Monge-Ampère partial differential equation[END_REF]. For a start, S M will be based on a finite difference scheme. Later on we will also propose a definition of S M based on a semi-Lagriangian scheme. Then, the filtered scheme S F is defined as

u n+1 j ≡ S F (u n ) j := S M (u n ) j + τ F S A (u n ) j -S M (u n ) j τ , (2.6)
where = τ,∆x > 0 is a parameter satisfying lim

(τ,∆x)→0 = 0. (2.7)
More hints on the choice of will be given later on.

The scheme is initialized in the standard way, i.e.

u 0 j := v 0 (x j ), ∀j ∈ Z. (2.8)
Now we make precise some requirements on S M , S A and the function F .

Definition of the monotone finite difference scheme S M Following Crandall and Lions [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF], we consider a finite difference scheme written as

u n+1 = S M (u n ) with S M (u n )(x) := u n (x) -τ h M (x, D -u n (x), D + u n (x)), (2.9)
with

D ± u(x) := ± u(x ± ∆x) -u(x) ∆x ,
where h M corresponds to a monotone numerical Hamiltonian that will be made precise below. We will denote also S M (u n ) j := S M (u n )(x j ). Therefore the scheme also reads, for all j ∈ Z, ∀n ≥ 0:

u n+1 j := u n j -τ h M (x j , D -u n j , D + u n j ), D ± u n j := ± u n j±1 -u n j ∆x . (2.10) (A3) Assumptions on S M
We will use the following assumptions throughout this paper:

(i) h M is a Lipschitz continuous function. (ii) (consistency) ∀x, ∀u, h M (x, u, u) = H(x, u). (iii) (monotonicity) for any functions u, v, u ≤ v =⇒ S M (u) ≤ S M (v).
In pratice condition (A3)-(iii) is only required at mesh points and the condition reads

u j ≤ v j , ∀j, ⇒ S M (u) j ≤ S M (v) j , ∀j.
(2.11)

At this stage, we notice that under condition (A3) the filtered scheme is " -monotone" in the sense that

u j ≤ v j , ∀j, ⇒ S M (u) j ≤ S M (v) j + τ, ∀j.
(2.12)

with → 0 as (τ, ∆x) → 0. This implies the convergence of the scheme by Barles-Souganidis convergence theorem (see [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF][START_REF] Abgrall | Construction of simple, stable and convergent hinge order scheme for steady first order hamiltonjacobi equation[END_REF]). Remark 2.1. Under assumption (i), the consistency property (ii) is equivalent to say that, for any v ∈ C 2 ([0, T ] × R), there exists a constant C M ≥ 0 independant of ∆x such that

h M (x, D -v(x), D + v(x)) -H(x, v x ) ≤ C M ∆x ∂ xx v ∞ .
(2.13)

The same statement holds true if (2.13) is replaced by the following consistency error estimate:

E S M (v)(t, x) := v(t + τ, x) -S M (v(t, .))(x) τ -v t (t, x) + H(x, v x (t, x))) ≤ C M τ ∂ tt v ∞ + ∆x ∂ xx v ∞ . (2.14)
Remark 2.2. Assuming (i), it is easily shown that the monotonicity property (iii) is equivalent to say that h M = h M (x, u -, u + ) satisfies, a.e. (x, u -, u + ) ∈ R 3 :

∂h M ∂u -≥ 0, ∂h M ∂u + ≤ 0, (2.15)
and the CFL condition

τ ∆x ∂h M ∂u -(x, u -, u + ) - ∂h M ∂u + (x, u -, u + ) ≤ 1. (2.16)
When using finite difference schemes, it is assumed that the CFL condition (2.16) is satisfied, and that can be written equivalently in the form

c 0 τ ∆x ≤ 1, (2.17)
where c 0 is a constant independant of τ and ∆x. Example 2.1. Let us consider the Lax-Friedrichs numerical Hamiltonian is

h M,LF (x, u -, u + ) := H(x, u -+ u + 2 ) - c 0 2 (u + -u -)
where c 0 > 0 is a constant. The scheme is consistant; it is furthermore monotone under the conditions max x,p |∂ p H(x, p)| ≤ c 0 , and c 0 τ ∆x ≤ 1. Definition of the high-order scheme S A : we consider an iterative scheme of "high-order" in the form u n+1 = S A (u n ), written as

S A (u n )(x) = u n (x) -τ h A (x, D k,-u n (x), . . . , D -u n (x), D + u n (x), . . . , D k,+ u n (x)),
where h A corresponds to a "high-order" numerical Hamiltonian, and D ,± u(x) := ± u(x± ∆x)-u(x) ∆x for = 1, . . . , k. To simplify the notations we may write (2.18) in the more compact form

S A (u n )(x) = u n (x) -τ h A x, D ± u n (x)) (2.18)
even if there is a dependency on in (D ,± u n (x)) =1,...,k .

(A4) Assumptions on S A :

We will use the following assumptions: (i) h A is a Lipschitz continuous function.

(ii) (high-order consistency) There exists k ≥ 2, for all ∈ [1, . . . , k], for any function v = v(t, x) of class C +1 , there exists C A, ≥ 0,

E S A (v)(t, x) := v(t + τ, x) -S A (v(t, .))(x) τ -v t (t, x) + H(x, v x (t, x))) (2.19) ≤ C A, τ ∂ +1 t v ∞ + ∆x ∂ +1 x v ∞ .
(2.20)

Here v x denotes the -th derivative of v w.r.t. x. Remark 2.3. The high-order consistency implies, for all ∈ [1, . . . , k], and for v ∈ C +1 (R),

h A (x, . . . , D -v, D + v, . . . ) -H(x, v x ) ≤ C A, ∂ +1 x v ∞ ∆x .
Example 2.2. (Centered scheme) A typical example with k = 2 is obtained with the centered TVD (Total Variation Diminishing) approximation in space and the Runge-Kutta 2nd order scheme in time (or Heun scheme):

S 0 (u n ) j := u n j -τ H(x j , u n j+1 -u n j-1 2∆x ), (2.21a)
and

S A (u) := 1 2 (u + S 0 (S 0 (u))). ( 2

.21b)

Of course there is no reason for the centered scheme to be stable (as it will be shown in the numerical section). Using a filter will help stabilize the scheme. A similar example with k = 3 can be obtained with any third order finite difference approximation in space and the TVD-RK3 scheme in time [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF].

Definition of the filter function F . We recall that Froese and Oberman's filter function used in [START_REF] Froese | Convergent filtered schemes for the Monge-Ampère partial differential equation[END_REF] is:

F (x) = sign(x) max(1 -||x| -1|, 0) =        x |x| ≤ 1. 0 |x| ≥ 2. -x + 2 1 ≤ x ≤ 2. -x -2 -2 ≤ x ≤ -1.
In the present work we define a new filter function simply as follows:

F (x) := x1 |x|≤1 = x if |x| ≤ 1, 0 otherwise. (2.22)
The idea of the present filter function is to keep the high-order scheme when

|h A -h M | ≤ (because then |S A -S M |/(τ ) ≤ 1 and S F = S M + τ F ( S A -S M τ
) ≡ S A ), whereas F = 0 and S F = S M if that bound is not satisfied, i.e., the scheme is simply given by the monotone scheme itself. Clearly the main difference is the discontinuity at x = -1, 1. 2.3. Convergence result. The following theorem gives several basic convergence results for the filtered scheme. Note that the high-order assumption (A4) will not be necessary to get the error estimates (i)-(ii). It will be only used to get a high-order consistency error estimate in the regular case (part (iii)). Globally the scheme will have just an O( √ ∆x) rate of convergence for just Lipschitz continuous solutions because the jumps in the gradient prevent high-order accuracy on the kinks.

Theorem 2.1. Assume (A1)-(A2), and let v 0 be bounded. We assume also that S M satisfies (A3), and |F | ≤ 1. Let u n denote the filtered scheme (2.6). Let v n j := v(t n , x j ) where v is the exact solution of (2.4). Assume

0 < ≤ c 0 √ ∆x (2.23)
for some constant c 0 > 0.

(i) The scheme u n satisfies the Crandall-Lions estimate

u n -v n ∞ ≤ C √ ∆x, ∀ n = 0, ..., N. (2.24)
for some constant C independent of ∆x.

(ii) (First order convergence for classical solutions.) If furthermore the exact solution v belongs to C 2 ([0, T ] × R), and ≤ c 0 ∆x (instead of (2.23)), then, we have

u n -v n ∞ ≤ C∆x, n = 0, ..., N, (2.25)
for some constant C independent of ∆x.

(iii) (Local high-order consistency.) Assume that S A is a high-order scheme satisfying (A4) for some k ≥ 2. Let 1 ≤ ≤ k and v be a C +1 function in a neighborhood of a point (t, x) ∈ (0, T )×R. Assume that

(C A,1 + C M ) v tt ∞ τ + v xx ∞ ∆x ≤ . (2.26)
Then, for sufficiently small t n -t, x j -x, τ , ∆x, it holds

S F (v n ) j = S A (v n ) j
and, in particular, a local high-order consistency error for the filtered scheme S F holds:

E S F (v n ) j ≡ E S A (v n ) j = O(∆x ) (the consistency error E S A is defined in (2.19)).
Proof. (i) Let w n+1 j = S M (w n ) j be defined with the monotone scheme only, with w 0 j = v 0 (x j ) = u 0 j . By definitions,

u n+1 j -w n+1 j = S M (u n ) j -S M (w n ) j + τ F .)
Hence, by using the monotonicity of S M , max

j |u n+1 j -w n+1 j | ≤ max j |u n j -w n j | + τ,
and by recursion, for n ≤ N , max

j |u n j -w n j | ≤ nτ ≤ T .
On the other hand, by Crandall and Lions [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF], an error estimate holds for the monotone scheme:

max j |w n j -v n j | ≤ C √ ∆x,
for some C ≥ 0. By summing up the previous bounds, we deduce max

j |u n j -v n j | ≤ C √ ∆x + T ,
and together with the assumption on , it gives the desired result.

(ii) Let E n j :=

v n+1 j -S M (v n ) j τ .
If the solution is C 2 regular with bounded second order derivatives, then the consistency error is bounded by

|E n j | ≤ C M (τ + ∆x). (2.27) Hence |u n+1 j -v n+1 j | = |S M (u n ) j -S M (v n ) j + τ E n j + τ F .)| ≤ u n -v n ∞ + τ E n ∞ + τ .
By recursion, for nτ ≤ T ,

u n -v n ∞ ≤ u 0 -v 0 ∞ + T ( max 0≤k≤N -1 E k ∞ + ).
Finally by using the assumption on , the bound (2.27) and the fact that τ = O(∆x) (using CFL condition (2.17)), we get the desired result.

(iii) To prove that S F (v n ) j = S A (v n ) j , one has to check that

|S A (v n ) j -S M (v n ) j | τ ≤ 1
as (τ, ∆x) → 0. By using the consistency error definitions,

|S A (v n ) j -S M (v n ) j | τ = v n+1 j -S M (v n ) j τ + v t (t n , x j ) + H(x j , v x (t n , x j )) - v n+1 j -S A (v n ) j τ + v t (t n , x j ) + H(x j , v x (t n , x j )) ≤ |E S M (v n ) j | + |E S A (v n ) j | ≤ (C A,1 + C M )(τ v tt ∞ + ∆x v xx ∞ )
Hence the desired result follows. Remark 2.4. Other approaches. It is already known from the original work of Osher and Shu [START_REF] Osher | High order essentially nonoscillatory schemes for hamilton-jacobi equations[END_REF] that it is possible to modify an ENO scheme in order to obtain a convergent scheme. For instance, if D ±,A u n j denotes a high-order finite difference derivative estimate (of ENO type), a projection on the first order finite difference derivative D ± u n j can be used, up to a controlled error (see in particular Remark 2.2 of [START_REF] Osher | High order essentially nonoscillatory schemes for hamilton-jacobi equations[END_REF]):

instead of D ±,A u n j , use P [D ± u n j ,M ∆x] (D ±,A u n j )
where P [a,b] (y) is the projection defined by:

P [a,b] (y) :=    y if a -b ≤ y ≤ a + b a -b if y ≤ a -b a + b if y ≥ a + b
and M > 0 is some constant greater than the expected value 1 2 |u xx (t n , x j )|. However, we emphasize that in our approach we do not consider a projection but a perturbation with a filter, which is sligthly different. Indeed, by using a projection into an interval of the form [a -M ∆x, a + M ∆x] where a = D ± u n i , numerical tests show that we may choose too often one of the extremal values a ± M ∆x which is then produces an overall too big error (worse than using the first order finite differences).

Following the present approach, we would rather advice to use, instead of D ±,A u n j , the value P[D ± u n j ,M ∆x] (D ±,A u n j )

where P[a,b] (y) is defined by: (We notice that (A2) is satisfied for hamiltonian functions such as (2.28).) Let [u] denote the P 1 -interpolation of u in dimension one on the mesh (x j ), i.e.

P[a,b] (y) := y if a -b ≤ y ≤ a + b a if y / ∈ [a -b, a + b]
x

∈ [x j , x j+1 ] ⇒ [u](x) := x j+1 -x ∆x u j + x -x j ∆x u j+1 . (2.30)
Then a monotone SL scheme can be defined as follows:

S M (u n ) j := min a∈A max b∈B [u n ] x j + τ f (x j , a, b) + τ (x j , a, b) .
(2.31)

A filtered scheme based on SL can then be defined by (2.6) and (2.31). Convergence result as well as error estimates could also be obtained in this framework. (For error estimates for the monotone SL scheme, we refer to [START_REF] Soravia | Estimates of convergence of fully discrete schemes for the Isaacs equation of pursuit-evasion differential games via maximum principle[END_REF][START_REF] Falcone | Semi-Lagrangian Approximation Schemes for Linear and Hamilton-Jacobi Equations[END_REF].)

2.4. Adding a limiter. The basic filter scheme (2.6) is designed to be of high-order where the solution is regular and when there is no viscosity aspects. However, for instance in the case of front propagation, it can be observed that the filter scheme may let small errors occur near extrema, when two possible directions of propagation occur in the same cell. This is the case for instance near a minima for an eikonal equation. In order to improve the scheme near extrema, we propose to introduce a limiter before doing the filtering process. Limiting correction will be needed only when there is some viscosity aspect (it is not needed for advection).

Let us consider the case of front propagation, i.e., equation of type (2.4), now with

H(x, v x ) = max a∈A f (x, a)v x (2.32)
(i.e., no distributive cost in the Hamiltonian function).

In the one-dimensional case, a viscosity aspect may occur at a minima detected at mesh point x i if min a f (x j , a) ≤ 0 and max a f (x j , a) ≥ 0.

(2.33)

In that case, the solution should not go below the local minima around this point, i.e., we want u n+1 j ≥ u min,j := min(u n j-1 , u n j , u n j+1 ), (2.34) and, in the same way, we want to impose that

u n+1 j ≤ u max,j := max(u n j-1 , u n j , u n j+1 ). ( 2

.35)

If we consider the high-order scheme to be of the form u n+1 j = u n j -τ h A (u n ), then the limiting process amounts to saying that

h A (u n ) j ≤ h max j := u n j -u min,j τ .
and

h A (u n ) j ≥ h min j := u n j -u max,j τ .
This amounts to define a limited hA such that, if (2.33) holds at mesh point x j , then hA (u n ) j := min max(h

A (u n ) j , h min j ), h max j .
and, otherwise,

hA j :≡ h A j .
Then the filtering process is the same, using hA instead of h A for the definition of the high-order scheme S F . For two dimensional equations a similar limiter could be developped in order to make the scheme more efficient at singular regions. However, for the numerical tests of the next section (in two dimensions) we will simply limit the scheme by using an equivalent of (2.34)-(2.35). Hence, instead of the scheme value u n+1 ij = S A (u n ) ij for the high-order scheme, we will update the value by

u n+1 ij = min(max(S A (u n ) ij , u min ij ), u max ij ), (2.36)
where

u min ij = min(u n ij , u n i±1,j , u n i,j±1
) and

u max ij = max(u n ij , u n i±1,j , u n i,j±1
). 2.5. How to choose the parameter : a simplified approach. The scheme should switch to high-order scheme when some regularity of the data is detected, and in that case we should have

S A (v) -S M (v) τ = h A (•) -h M (•) ≤ 1.
In a region where a function v = v(x) is regular enough, by using Taylor expansions, zero order terms in h A (x, D ± v) and h M (x, D ± v) vanish (they are both equal to H(x, v x (x))) and it remains an estimate of order ∆x. More precisely, by using the high-order property (A4) we have

h A (x j , D ± v j ) = H(x j , v x (x j )) + O(∆x 2 ).
On the other hand, by using Taylor expansions,

Dv ± j = v x (x j ) ± 1 2 v xx (x j )∆x + O(∆x 2 ),
Hence, denoting h M = h M (x, u -, u + ), it holds at points where h M is regular,

h M (x j , Dv - j , Dv + j ) = H(x j , v x (x j )) + 1 2 v xx (x j ) ∂h M j ∂u + - ∂h M j ∂u -+ O(∆x 2 ).
Therefore,

|h A (v) -h M (v)| = 1 2 |v xx (x j )| ∂h M j ∂u + - ∂h M j ∂u -∆x + O(∆x 2 ).
Hence we will make the choice to take roughly such that

1 2 |v xx (x j )| ∂h M j ∂u + - ∂h M j ∂u -∆x ≤ (2.

37)

(where h M j = h M (x j , v x (x j ), v x (x j ))). Therefore, if at some point x j (2.37) holds, then the scheme will switch to the high-order scheme. Otherwise, when the expectations from h M and h A are different enough, the scheme will switch to the monotone scheme.

In conclusion we have upper and lower bound for the switching parameter :

• Choose ≤ c 0 √ ∆x for some constant c 0 > 0 in order that the convergence and error estimate result holds (see Theorem 2.1).

• Choose ≥ c 1 ∆x, where c 1 is sufficiently large. This constant should be choosen roughly such that

1 2 v xx ∞ ∂h M ∂u + (., v x , v x ) - ∂h M ∂u -(., v x , v x ) ∞ ≤ c 1 .
where the range of values of v x and v xx can be estimated, in general, from the values of (v 0 ) x , (v 0 ) xx and the Hamiltonian function H. Then the scheme is expected to switch to the high-order scheme where the solution is regular.

Numerical tests.

In this section we present several numerical tests in one and two dimensions. Unless otherwise indicated, the filtered scheme will refer to the scheme where the high-order Hamiltonian is the centered scheme in space (see Remark 2.2), with Heun (RK2) scheme discretization in time (see in particular Eqs. (2.21a)-(2.21b)). Hereafter this scheme will be referred as the "centered scheme".

The monotone finite difference scheme and function h M will be made precise for each example.

For the filtered scheme, unless otherwise precised, the switching coefficient = 5∆x. will be used. In practice we have numerically observed that taking = c 1 ∆x with c 1 sufficiently large does not much change the numerical results in the following tests. All the tested filtered schemes (apart from the steady and obstacle equations) enters in the convergence framework of the previous section, so in particular there is a theoretical convergence of order √ ∆x under the usual CFL condition.

In the tests, the filtered scheme will be in general compared to a second order ENO scheme (for precise definition, see Appendix A), as well as the centered (a priori unstable) scheme without filtering.

In several cases, local error in the L 2 norms are computed in some subdomain D, which, at a given time t n , corresponds to

e L 2 loc := {i, xi∈D} |v(t n , x i ) -u n i | 2 1/2
The first two numerical examples deal with one-dimensional HJ equations, examples 3 and 4 are concerned with two-dimensional HJ equations, and the last three examples will concern a one-dimensional steady equation and two nonlinear one-dimensional obstacle problems.

Example 1. Eikonal equation. We consider the case of

v t + |v x | = 0, t ∈ (0, T ), x ∈ (-2, 2), (3.1) v(0, x) = v 0 (x) := max(0, 1 -x 2 ) 4 , x ∈ (-2, 2). (3.2)
In Table 3.1, we compare the filtered scheme (with = 5∆x) with the centered scheme and the ENO second order scheme, with CFL = 0.37 and terminal time T = 0.3. For the filtered scheme, the monotone hamiltonian used is h

M (x, v -, v + ) := max(v -, -v + ).
As expected, we observe that the centered scheme alone is unstable. On the other hand, the filtered and ENO schemes are numerically comparable in that case, and second order convergent (the results are similar for the L 1 and the L ∞ errors).

Then, in Table 3.2, we consider the same PDE but with the following reversed initial data:

ṽ0 (x) := -max(0, 1 -x 2 ) 4 , x ∈ (-2, 2). (3.3)
In that case the centered scheme alone is unbounded. The filtered scheme (with = 5∆x) is second order. However, here, the limiter correction as described in section (2.4) get second order behavior. We also observe that the filtered scheme gives better results than the ENO scheme. (We have also numerically tested the ENO scheme with the same limiter correction but it does not improve the behavior of the ENO scheme alone).

In conclusion, this first example shows firstly, that the filtered scheme can stabilize an otherwise unstable scheme, and secondly that it can give the desired second order behavior. Example 2. Burger's equation.

In this example an HJ equivalent of the nonlinear Burger's equation is considered: with Dirichlet boundary condition on (-2, 2). Exact solution is known. 1 . In order to test highorder convergence we have considered the smoother initial data which is the one obtained from (3.4) at time t 0 := 0.1, i.e. :

v t + 1 2 |v x | 2 = 0, t > 0, x ∈ (-2, 2) (3.4a) v(0, x) = v 0 (x) := max(0, 1 -x 2 ), x ∈ (-2 , 2) (3.4b)
w t + 1 2 |w x | 2 = 0, t > 0, x ∈ (-2, 2). (3.5a) w(0, x) := v(t 0 , x), x ∈ (-2, 2), (3.5b)
with exact solution w(t, x) = v(t + t 0 , x).

An illustration is given in Fig. 3.3. For the filtered scheme, the monotone hamiltonian used is

h M (x, v -, v + ) := 1 2 (v -) 2 1 v ->0 + 1 2 (v + ) 2 1 v + <0 . Errors are given in Table (3.
3), using CFL=0.37 and terminal time T = 0.3. In conclusion we observe numerically that the filtered scheme keeps the good behavior of the centered scheme (here stable and almost second order).

Example 3. 2D rotation. We now apply filtered scheme to an advection equation in two dimensions:

1 It holds v(t, x) = (max(0, 1 -|x|)) 2 2t 1 {t> 1 2 } + (1 -2t) 2 -|x| 2 1 -2t 1 {1≥|x|≥1-2t} .
v t -yv x + xv y = 0, (x, y) ∈ Ω, t > 0, (3.6) v(0, x, y)
= v 0 (x, y) := 0.5 -0.5 max(0, 1 -(x -1) 2 -y 2 1 -r 2 0 ) 4 (3.7)
where Ω := (-A, A) 2 (with A = 2.5), r 0 = 0.5 and with Dirichlet boundary condition v(t, x) = 0.5, x ∈ ∂Ω. This initial condition is regular and such that the level set v 0 (x, y) = 0 corresponds to a circle centered at (1, 0) and of radius r 0 .

In this example the monotone numerical Hamiltonian is defined by

h M (u - x , u + x , u - y , u + y ) := max(0, f 1 (a, x, y))u - x + min(0, f 1 (a, x, y))u + x (3.8)
+ max(0, f 2 (a, x, y))u - y + min(0, f 2 (a, x, y))u + y and the high-order scheme is the centered finite difference scheme in both spacial variables, and RK2 in time. The filtered scheme is otherwise the same as (2.6). However it is necessary to use a greater constant c 1 is the choice = c 1 ∆x in order to get (global) second order errors. We have used here = 20∆x. On the other hand the CFL condition is

µ := c 0 ( τ ∆x + τ ∆y ) ≤ 1, (3.9)
where here c 0 = 2.5 (an upper bound for the dynamics in the considered domain Ω). In this test the CFL number is µ := 0.37.

Results are shown in Table 3.4 for terminal time time T := π/2. Although the centered scheme is a priori unstable, in this example it is numerically stable and of second order. So we observe that the filtered scheme keep this good behavior and is also or second order (ENO scheme gives comparable results here). The zero-level set of v 0 corresponds to two separates circles or radius r 0 and centered in A = (1, 0) and B = (-1, 0) respectively. Dirchlet boundary conditions are used as the previous example.

The monotone hamiltonian h M used in the filtered scheme is in Lax-Friedrichs form:

h M (x, u - 1 , u + 1 , u - 2 , u + 2 ) = H(x, u - 1 + u + 1 2 , u - 2 + u + 2 2 ) - C x 2 (u + 1 -u - 1 ) - C y 2 (u + 2 -u - 2 ), (3.12)
where, here, C x = C y = 1. We used the CFL condition µ = 0.37 as in (3.9). Also, the simple limiter (2.36) was used for the filtered scheme as described in Section 2.4. It is needed in order to get a good second order behavior at extrema of the numerical solution. The filter coefficient is set to = 20∆x as in the previous example. Numerical results are given in Table 3.5, showing the global L 2 errors for the filtered scheme, the centered scheme, and a second order ENO scheme, at time t = 0.6. We observe that the centered scheme has some unstabilities for fine mesh, while the filtered performs as expected.

Example 5 Steady eikonal equation. We consider a steady eikonal equation with Dirichlet boundary condition, which is taken from Abgrall [START_REF] Abgrall | Construction of simple, stable and convergent hinge order scheme for steady first order hamiltonjacobi equation[END_REF]:

|v x | = f (x) x ∈ (0, 1), (3.13a) v(0) = v(1) = 0, (3.13b)
where f (x) = 3x 2 + a, with a = The steady solution for (3.13) can be considered as the limit lim t→∞ v(t, x) where v is the solution of the time marching corresponding form:

1-2x 3 0 2x0-1 and x 0 = 3 √ 2+2 4 3 √ 2 . Exact solution is known: v(x) := x 3 + ax x ∈ [0, x 0 ], 1 + a -ax -x 3 x ∈ [x 0 , 1]. ( 3
v t + |v x | = f (x) x ∈ (0, 1), t > 0, (3.15a) v(t, 0) = v(t, 1) = 0, t > 0. (3.15b)
In this example the upwind monotone scheme is used:

h M (.) j := u n+1 j -u n j τ -max{ u n j -u n j-1 ∆x , u n j -u n j+1 ∆x } -τ f (x j ),
the high-order scheme will be the centered scheme, and the filtered scheme (2.6) will be used with = 5∆x. The iterations are stopped when the difference between too successive time step is small enough or a fixed number of iterations is passed, i.e., in this example,

u n+1 -u n L ∞ := max i |u n+1 i -u n i | ≤ 10 -6 or n ≥ N max := 5000. (3.16)
As analyzed in [START_REF] Bokanowski | Value iteration convergence of -monotone schemes for stationary hamilton-jacobi equations[END_REF] for -monotone schemes, for a given mesh step, even if the iterations may not converge (because of the non monotony of the scheme), it can be shown to be close to a fixed point after enough iterations. Example 6 Advection with an obstacle. Here we consider an obstacle problem, which is taken from [START_REF] Bokanowski | A discontinuous galerkin scheme forefront propagation with obstacle[END_REF]:

min(v t + v x , v -g(x)) = 0, t > 0, x ∈ [-1, 1], (3.17) v 0 (x) = 0.5 + sin(πx) x ∈ [-1, 1], (3.18)
together with periodic boundary condition. The obstacle function is g(x) := sin(πx). In this case exact solution is given by:

v(t, x) :=    max(v 0 (x -at), g(x)) if t < 1 3 max(v 0 (x -at), g(x), -1 x∈ [0.5,1] ) if t ∈ [ 1 3 , 5 6 ], max(v 0 (x -at), g(x), 1 x∈ [-1,t-5 6 ]∪[0.5,1] ) if t ∈ [ 5 6 , 1], (3.19)
Results are given in with periodic boundary condition on (-1, 1) and g(x) = sin(πx). In this case the exact solution is given by: v(t, x) = max(v(t, x), g(x)). Results are given in Table 3.8 for terminal time T = 0.2. Plots are also shown in Figure 3.8 for different times (for t ≥ 1 3 solution remains unchanged). 4. Conclusion. We propose a "filtered" scheme which behaves as a high-order scheme when the solution is smooth and as a low order monotone scheme otherwise. It has a simple presentation that is easy to implement. Rigorous error bounds hold, of the same order as the Crandall-Lions estimates in √ ∆x where ∆x is the mesh size. In the case the solution is smooth a high-order consistency error estimate also holds. We show on several numerical examples the ability of the scheme to stabilize an otherwise unstable scheme, and also we observe a precision similar to a second order ENO scheme on basic linear and non linear examples. On going works concern the application of the present approach to some front propagation equations.

2. 1 .

 1 Setting of the problem. Let us denote by | • | the Euclidean norm on R d (d ≥ 1). The following classical assumptions will be considered in the sequel of this paper: (A1) v 0 is Lipschitz continuous function i.e. there exist L 0 > 0 such that for every x , y ∈ R d , |v 0 (x) -v 0 (y)| ≤ L 0 |x -y|.

Fig. 2 . 1 .

 21 Fig. 2.1. Froese and Oberman's filter (left), new filter (right)

Remark 2 . 5 .

 25 Filtered semi-Lagrangian scheme. Let us consider the case of H(x, p) := min b∈B max a∈A {-f (x, a, b).p -(x, a, b)}, (2.28) where A ⊂ R m and B ⊂ R n are non-empty compact sets (with m, n ≥ 1), f : R d × A × B → R d and : R d × A × B → R are Lipschitz continuous w.r.t. x: ∃L ≥ 0, ∀(a, b) ∈ A × B, ∀x, y: max(|f (x, a, b) -f (y, a, b)|, | (x, a, b) -(y, a, b)|) ≤ L|x -y|. (2.29)

Fig. 3 . 1 .

 31 Fig. 3.1. (Example 1) With initial data (3.2) (left), and plots at time T = 0.3 with centered scheme -middle -and filtered scheme -right, using M = 160 mesh points.

Fig. 3 . 2 .

 32 Fig. 3.2. (Example 1) With initial data (3.3) (left), and plots at time T = 0.3 with centered scheme -middle -and filtered scheme -right, using M = 160 mesh points.

13

 13

Fig. 3 . 3 .

 33 Fig. 3.3. (Example 2) Plots at t = 0 and t = 0.3 with the filtered scheme.

Example 4 .-r 2 0) 4 ,-r 2 0 ) 4 .

 42424 Eikonal equation. In this example we consider the eikonal equationv t + |∇v| = 0, (x, y) ∈ Ω, t > 0 (3.10) in the domain Ω := (-3, 3) 2 . The initial data is defined by v0(x, y) = (3.11) 0.5 -0.5 max max(0, 1 -(x -1) 2 -y 2 1 max(0, 1 -(x + 1) 2 -y 2 1 filtered centered ENO M x N y L 2 error order L 2 error order L 2 error order 20 77E-02 1.98 3.77E-02 1.98 2.04E-01 1.19 160 160 9.40E-03 2.00 9.40E-03 2.00 5.50E-02 1.89 320 320 2.34E-03 2.01 2.34E-03 2.01 1.29E-02 2.10

Table 3. 4 (

 4 Example 3) Global L 2 errors for the filtered scheme, centered and second order ENO schemes (with CFL 0.37).

Fig. 3 . 4 .

 34 Fig. 3.4. (Example 3) Filtered scheme, plots at time t = 0 (left) and t = π/2 (rigth) with M = 80 mesh points.

Fig. 3 . 5 .

 35 Fig. 3.5. (Example 4) Plots at times t = 0 (top) and t = π/2 (bottom) for the filtered scheme with M = 50 mesh points. The figures to the right represent the 0-level sets.

Fig. 3 . 6 .

 36 Fig. 3.6. (Example 5) Filtered scheme for a steady equation, with M = 50 mesh points.

Fig. 3 . 7 . 5 . 7

 3757 Fig. 3.7. (Example 6) Plots at T=0(initial data), T=0.3, T=0.5.

  (3.22) where v is the solution of the Eikonal equationv t +|v x | = 0. The formula v(t, x) = min y∈[x-t,x+t] v 0 (y) holds, which simplifies to v(t, x) :=    v 0 (x + t) if x < -0.5 -t -0.5 if x ∈ [-0.5 -t, -0.5 + t], min(v 0 (x -t), v 0 (x + t)) if x ≥ -0.5 + t,(3.23)

Fig. 3 . 8 .

 38 Fig. 3.8. (Example 7) Plots at times t = 0, t = 0.2 and t = 0.4. The dark line is the numerical solution, similar to the exact solution, and the ligth line is the obstacle function.

Table 3

 3

	.1

Table 3 .

 3 2(Example 1 with initial data (3.3).) L 2 errors for filtered scheme, centered scheme, and ENO second order scheme.

Table 3 .

 3 3(Example 2) L 2 errors for filtered scheme, centered scheme, and ENO second order scheme.

Table 3 .

 3 5(Example 4) Global L 2 errors for filtered scheme, centered and second order ENO schemes.

Table 3 .

 3 6(Example 5) Global errors for filtered scheme, compared with the centered (unstable) scheme, and a filtered ENO scheme.

		filtered				centered		filtered + ENO
	M	error	order error order		error	order
	50 2.16E-03		-		NaN	-			5.29E-03	-
	100 7.14E-04 1.60		NaN	-			1.35E-03 1.97
	200 2.17E-04 1.72		NaN	-			3.42E-04 1.98
	400 6.32E-05 1.78		NaN	-			8.61E-05 1.99
	800 2.17E-05 1.54		NaN	-			2.16E-05 2.00
		1.5					Exact(red), Scheme(blue)
		1
		0.5
		0
		-0.5	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Table 3 .

 3 [START_REF] Crandall | Two approximations of solutions of Hamilton-Jacobi equations[END_REF], for terminal time t = 0.5. Errors are computed away from singular points, i.e., in the region [-1, 1] \ ∪ i=1,3 [s i -δ, s i + δ] (where s 1 = -0.1349733, s 2 = 0.5 and s 3 = 2/3 are the three singular points. Filtered scheme is numerically of second order (ENO gives comparable results here).

	Errors	filtered = 5∆x	centered	ENO2
	M	N	error	order	error	order	error	order
	40	20 7.93E-03 2.03 1.63E-02 1.54 2.14E-02 1.59
	80	40 1.84E-03 2.10 2.98E-02 -0.87 7.75E-03 1.46
	160 80 3.92E-04 2.24 1.46E-02 1.03 1.07E-03 2.86
	320 160 9.67E-05 2.02 8.02E-03 0.86 2.72E-04 1.97
	640 320 2.40E-05 2.01 4.10E-03 0.97 6.92E-05 1.98

Table 3 .

 3 7(Example 6) L ∞ errors away from singular points, for filtered scheme, centered scheme, and second order ENO scheme.

Table 3 . 8

 38 Filtered scheme and ENO scheme at time t = 0.2

											Errors	filtered	ENO2
											M	error	order	error	order
											40	3.74E-03	-	6.85E-03	-
											80	6.26E-04 2.58 2.12E-03 1.69
											160	1.13E-04 2.47 6.80E-04 1.64
											320	2.26E-05 2.32 2.18E-04 1.64
											640	5.50E-06 2.04 6.96E-05 1.65
	-1 -1.5	-0.8	-0.6	-0.4	-0.2	0	0.2	0.4	0.6	0.8	1

The authors wish to acknowledge the support obtained by the following grants: AFOSR Grant no. FA9550-10-1-0029, and by the EU under the 7th Framework Programme Marie Curie Initial Training Network "FP7-PEOPLE-2010-ITN", SADCO project, GA number 264735-SADCO

Appendix A. An essentially non-oscillatory (ENO) scheme of second order. We recall here a simple ENO method of order two based on the work of Osher and Shu [START_REF] Osher | High order essentially nonoscillatory schemes for hamilton-jacobi equations[END_REF] for Hamilton Jacobi equation (the ENO method was designed by Harten et al. [START_REF] Harten | Uniformly high order essentially non-oscillatory schemes[END_REF] for the approximation solution of non-linear conservation law).

Let m be the minmod function defined by

(other functions can be considered such as m(a, b) = a if |a| ≤ |b| and m(a, b) = b otherwise). Let D ± u j = ±(u j±1 -u j )/∆x and

Then the right and left ENO approximation of the derivative can be defined by

and the ENO (Euler forward) scheme by S 0 (u) j := u j -τ h M (x j , Du j , D+ u j ).

The corresponding RK2 scheme can then be defined by S(u) = 1 2 (u + S 0 (S 0 (u))).


Pierre-Alain Fouque
email: pierre-alain.fouque@ens.fr

Nicolas Guillermin
email: nicolas.guillermin@m4x.org

Delphine Leresteux
email: delphine.leresteux@dga.defense.gouv.fr

Mehdi Tibouchi
email: tibouchi.mehdi@lab.ntt.co.jp

Jean-Christophe Zapalowicz
email: jean-christophe.zapalowicz@inria.fr

Attacking RSA-CRT Signatures with Faults on Montgomery Multiplication

Keywords: Fault Attacks, Montgomery Multiplication, RSA-CRT, PSS

In this paper, we present several efficient fault attacks against implementations of RSA-CRT signatures that use modular exponentiation algorithms based on Montgomery multiplication. They apply to any padding function, including randomized paddings, and as such are the first fault attacks effective against RSA-PSS.

The new attacks work provided that a small register can be forced to either zero, or a constant value, or a value with zero high-order bits. We show that these models are quite realistic, as such faults can be achieved against many proposed hardware designs for RSA signatures.

Introduction

The RSA signature scheme is one of the most used schemes nowadays. An RSA signature is computed by applying some encoding function to the message, and raising the result to d-th power modulo N , where d and N are the private exponent and the public modulus respectively. This modular exponentiation is the costlier part of signature generation, so it is important to implement it efficiently. A very commonly used speedup is the RSA-CRT signature generation, where the exponentiation is carried out separately modulo the two factors of N , and the results are then recombined using the Chinese Remainder Theorem. However, when unprotected, RSA-CRT signatures are vulnerable to the so-called Bellcore attack first introduced by Boneh et al. in [START_REF] Boneh | On the importance of checking cryptographic protocols for faults[END_REF], and later refined in multiple publications such as [START_REF] Yen | Hardware fault attack on RSA with CRT revisited[END_REF]: an attacker who knows the padded message and is able to inject a fault in one of the two half-exponentiations can factor the public modulus using a faulty signature with a simple GCD computation.

Many workarounds have been proposed to patch this vulnerability, including extra computations and sanity checks of intermediate and final results. A recent taxonomy of these countermeasures is given in [START_REF] Rivain | Securing RSA against fault analysis by double addition chain exponentiation[END_REF]. The simplest countermeasure may be to verify the signature before releasing it. This is reasonably cheap if the public exponent e is small and available in the signing device. In some cases, however, e is not small, or even not given-e.g. the JavaCard API does not provide it [22]. Another approach is to use an extended modulus. Shamir's trick [START_REF] Shamir | Improved method and apparatus for protecting public key schemes from timing and fault attacks[END_REF] was the first such technique to be proposed; later refinements were suggested that also protect CRT recombination when it is computed using Garner's formula [START_REF] Blömer | A new CRT-RSA algorithm secure against Bellcore attacks[END_REF][START_REF] Ciet | Practical fault countermeasures for Chinese remaindering based cryptosystems[END_REF][START_REF] Vigilant | RSA with CRT: A new cost-effective solution to thwart fault attacks[END_REF][START_REF] Coron | Fault attacks and countermeasures on Vigilant's RSA-CRT algorithm[END_REF]. Finally, yet another way to protect RSA-CRT signatures against faults is to use redundant exponentiation algorithms, such as the Montgomery Ladder. Papers including [START_REF] Giraud | An RSA implementation resistant to fault attacks and to simple power analysis[END_REF][START_REF] Rivain | Securing RSA against fault analysis by double addition chain exponentiation[END_REF] propose such countermeasures. Regardless of the approach, RSA-CRT fault countermeasures tend to be rather costly: for example, Rivain's countermeasure [START_REF] Rivain | Securing RSA against fault analysis by double addition chain exponentiation[END_REF] has a stated overhead of 10% compared to an unprotected implementation, and is purportedly more efficient than previous works including [START_REF] Giraud | An RSA implementation resistant to fault attacks and to simple power analysis[END_REF][START_REF] Vigilant | RSA with CRT: A new cost-effective solution to thwart fault attacks[END_REF].

Relatedly, while Boneh et al.'s original fault attack does not apply to RSA signatures with probabilistic encoding functions, some extensions of it were proposed to attack randomized ad hoc padding schemes such as ISO 9796-2 and EMV [START_REF] Coron | Fault attacks on RSA signatures with partially unknown messages[END_REF][START_REF] Coron | Fault attacks against EMV signatures[END_REF]. However, Coron and Mandal [START_REF] Coron | PSS is secure against random fault attacks[END_REF] were able to prove that Bellare and Rogaway's padding scheme RSA-PSS [START_REF] Bellare | Probabilistic signature scheme[END_REF] is secure against random faults in the random oracle model. In other words, if injecting a fault on the half-exponentiation modulo the second factor q of N produces a result that can be modeled as uniformly distributed modulo q, then the result of such a fault cannot be used to break RSA-PSS signatures. It is tempting to conclude that using RSA-PSS should enable signers to dispense with costly RSA-CRT countermeasures. In this paper, we argue that this is not necessarily the case.

Our contributions. The RSA-CRT implementations targeted in this paper use the state-of-the-art modular multiplication algorithm due to Montgomery [START_REF] Montgomery | Modular multiplication without trial division[END_REF], which avoids the need to compute actual divisions on large integers, replacing them with only multiplications and bit shifts. A typical implementation of the Montgomery multiplication algorithm will use small registers to store precomputed values or short integer variables throughout the computation. The size of these registers varies with the architecture, from a single bit in certain hardware implementations to 16 bits, 32 bits or more in software. This paper presents several fault attacks on these small registers during Montgomery multiplication, that cause the result of one of the half-exponentiations to be unusually small. The factorization of N can then be recovered using a GCD, or an approximate common divisor algorithm such as [START_REF] Howgrave-Graham | Approximate integer common divisors[END_REF][START_REF] Chen | Faster algorithms for approximate common divisors[END_REF][START_REF] Cohn | Approximate common divisors via lattices[END_REF].

We consider three models of faults on the small registers. In the first model, one register can be forced to zero. In that case, we show that causing such a fault in the inverse Montgomery transformation of the result of a half-exponentiation, or a few earlier consecutive Montgomery multiplications, yields a faulty signature which is a multiple of the corresponding factor q of N . Hence, we can factor N by taking a simple GCD. In the second model, another register can be forced to some (possibly unknown) constant value throughout the inverse Montgomery transformation of the result of a half-exponentiation, or a few earlier consecutive Montgomery multiplications. A faulty signature in this model is a close multiple of the corresponding factor q of N , and we can thus factor N using an approximate common divisor algorithm. Finally, the third model makes it possible to force some of the higher-order bits of one register to zero. We show that, while injecting one such fault at the end of the inverse Montgomery transformation results in a faulty signature that isn't usually close enough to a multiple of q to reveal the factorization of N on its own, a moderate number of faulty signatures (a dozen or so) obtained using that process are enough to factor N .

The RSA padding scheme used for signing, whether deterministic or probabilistic, is irrelevant in our attacks. In particular, RSA-PSS implementations are also vulnerable. Of course, this does not contradict the security result due to Coron and Mandal [START_REF] Coron | PSS is secure against random fault attacks[END_REF], as the faults we consider are strongly non-random. Our results do suggest, however, that exponentiation algorithms based on Montgomery multiplication are quite sensitive to a very realistic type of fault attacks and that using RSA-CRT countermeasures is advisable even for RSA-PSS.

Organization of the paper. In §2, we recall some background material on the Montgomery multiplication algorithm, on modular exponentiation techniques, and on RSA-CRT signatures. Our new attacks are then described in § §3-5, corresponding to three different fault models: null faults, constant faults, and zero high-order bits faults. Finally, in §6, we discuss the applicability of our fault models to concrete hardware implementations of RSA-CRT signatures, and find that many proposed designs are vulnerable.

Preliminaries

Montgomery multiplication

Proposed by Montgomery in [START_REF] Montgomery | Modular multiplication without trial division[END_REF], the Montgomery multiplication algorithm provides a fast way method for computing modular multiplications and squarings. Indeed, the Montgomery multiplication algorithm only uses multiplications, additions and shifts, and its cost is about twice that of a simple multiplication (compared to 2.5 times for a multiplication and a Barett reduction), without imposing any constraint on the modulus.

Usually, one of two different techniques is used to compute Montgomery multiplication: either Separate Operand Scanning (SOS), or Coarsely Integrated Operand Scanning (CIOS). Consider a device whose processor or coprocessor architecture has r-bit registers (typically r = 1, 8, 16, 32 or 64 bits). Let b = 2 r , q be the (odd) modulus with respect to which multiplications are carried out, k the number of r-bit registers used to store q, and R = b k , so that q < R and gcd(q, R) = 1. The SOS variant consists in using the Montgomery reduction after the multiplication: for an input A such that A < Rq, it computes Mgt(A) ≡ AR -1 (mod q), with 0 ≤ Mgt(A) < q. The CIOS mixes the reduction algorithm with the previous multiplication step: considering x and y with xy < Rq, it computes CIOS(x, y) = xyR -1 mod q with CIOS(x, y) < q.

Figure 2 presents the main steps of the CIOS variant, which will be used thereafter. However, replacing the CIOS by the SOS or any other variant proposed in [START_REF] Koç | Analyzing and comparing Montgomery multiplication algorithms[END_REF] does not protect against any of our attacks.

Exponentiation algorithms using Montgomery multiplication

Montgomery reduction is especially interesting when used as part of a modular exponentiation algorithm. A large number of such exponentiation algorithms are known, including the Square-and-Multiply algorithm from either the least or the most significant bit of the exponent, the Montgomery Ladder (used as a side-channel countermeasure against cache analysis, branch analysis, timing analysis and power analysis), the Squareand-Multiply k-ary algorithm (which boasts greater efficiency thanks to fewer multiplications), etc. The first three exponentiation algorithms will be considered in this paper, and two of those are detailed in Figure 3. Note that using the Montgomery multiplications inside any exponentiation algorithm requires all variables to be in Montgomery representation (x = xR mod q is the Montgomery representation of x) before applying

1: function SignRSA-CRT(m) 2: M ← µ(m) ∈ ZN ⊲ message encoding 3: Mp ← M mod p 4: Mq ← M mod q 5: Sp ← M dp p mod p 6: Sq ← M dq q mod q 7: t ← Sp -Sq 8: if t < 0 then t ← t + p 9: S ← Sq + (t • π) mod p • q 10:
return S Fig. 1. RSA-CRT signature generation with Garner's recombination. The reductions dp, dq modulo p -1, q -1 of the private exponent are precomputed, as is π = q -1 mod p.

1: function CIOS(x, y) 2: a ← 0 3: y0 ← y mod b 4: for j = 0 to k -1 do 5: a0 ← a mod b 6: uj ← (a0 + xj • y0) • q ′ mod b 7: a ← a + xj • y + uj • q b 8: if a ≥ q then a ← a -q 9:
return a

; q ′ = -q -1 mod b is precomputed. The re- turned value is (xy • b -k mod q). Since b = 2 r , the division is a bit shift.
the exponentiation process. In line 2 of each algorithm from Figure 3, the message is transformed into Montgomery representation by computing CIOS(x, R 2 ) = xR 2 R -1 mod q = x. At the end, the very last CIOS call allows to revert to the classical representation by performing a Montgomery reduction:

CIOS( Ā, 1) = ( Ā•1)R -1 mod q = ARR -1 mod q = A.
Finally the other CIOS steps compute the product in Montgomery representation: CIOS( Ā, B) = (AR)(BR)R -1 mod q = AB.

RSA-CRT signature generation

Let N = pq be a n-bit RSA modulus. The public key is denoted by (N, e) and the associated private key by (p, q, d). For a message M to be signed, we note S = m d mod N the corresponding signature, where m is deduced from M by an encoding function, possibly randomized. A well-known optimization of this operation is the RSA-CRT which takes advantage of the decomposition of N into prime factors. By replacing a full exponentiation of size n by two n/2, it divides the computational cost by a factor of around 4. Therefore RSA-CRT is almost always employed: for example, OpenSSL as well as the JavaCard API [22] use it.

Recovering S from its reductions S p and S q modulo p and q can be done either by the usual CRT reconstruction formula (1) below, or using the recombination technique (2) due to Garner:

S = (S q • p -1 mod q) • p + (S p • q -1 mod p) • q mod N.
(1)

S = S q + q • (q -1 • (S p -S q ) mod p). (2)
Square-and-Multiply LSB Montgomery Ladder 1: function ExpLSB(x, e, q) 2:

x ← CIOS(x, R 2 mod q) 3:

A ← R mod q 4: for i = 0 to t do 5:

if ei = 1 then 6:

A ← CIOS(A, x) 7:

x ← CIOS(x, x) 8:

A ← CIOS(A, 1) 9:

return A 1: function ExpLadder(x, e, q) 2:

x ← CIOS(x, R 2 mod q) 3:

A ← R mod q 4:

for i = t down to 0 do 5:

if ei = 0 then 6:

x ← CIOS(A, x) 7:

A ← CIOS(A, A) 8:

else if ei = 1 then 9:

A ← CIOS(A, x) 10:

x ← CIOS(x, x) 11:

A ← CIOS(A, 1) 12:

return A Fig. 3. Two of the exponentiation algorithms considered in this paper. In each case, e0, . . . , et are the bits of the exponent e (from the least to the most significant), b is the base in which computations are carried out (gcd(b, q) = 1) and R = b k .

Garner's formula (2) does not require a reduction modulo N , which is interesting for efficiency reasons and also because it prevents certain fault attacks [START_REF] Brier | Modulus fault attacks against RSA-CRT signatures[END_REF]. On the other hand, it does require an inverse Montgomery transformation S q = CIOS( Sq , 1), whereas that step is not necessary for formula [START_REF] Bellare | Probabilistic signature scheme[END_REF], as it can be mixed with the multiplication with q -1 mod p. This is an important point, as some of our attacks specifically target the inverse Montgomery transformation. The main steps of the RSA-CRT signature generation with Garner's recombination are recalled in Figure 1.

Null Faults

We first consider a fault model in which the attacker can force the register containing the precomputed value q ′ = (-q mod b) to zero in certain calls to the CIOS algorithm during the computation of S q .

Under suitable conditions, we will see that such faults can cause the qpart of the signature to be erroneously evaluated as S q = 0, which makes it possible to retrieve the factor q of N from one such faulty signature S, as q = gcd( S, N ).

Attacking CIOS(A, 1)

Suppose first that the fault attacker can force q ′ to zero in the very last CIOS computation during the evaluation of S q , namely the computation of CIOS(A, 1). In that case, the situation is quite simple. Theorem 1. A faulty signature S generated in this fault model is a multiple of q (for any of the exponentiation algorithms considered herein and regardless of the encoding function involved, probabilistic or not).

Proof. The faulty value q ′ = 0 causes all of the variables u in the CIOS loop to vanish; indeed, for j = 0, . . . , k -1, they evaluate to:

u j = (a 0 + A j • 1) • q ′ mod 2 r = 0.
As a result, the value S q computed by this CIOS loop can be written as:

S q = • • • A 0 • 2 -r + A 1 • 2 -r + • • • + A k-1 • 2 -r .
Now, the values A j are r-words, i.e. 0 ≤ A j ≤ 2 r -1. It follows that each of the integer divisions by 2 r evaluate to zero, and hence S q = 0. As a result, the faulty signature S is a multiple of q as stated.

⊓ ⊔

It is thus easy to factor N with a single faulty signature S, by computing gcd( S, N ). Note also that if this last CIOS step is computed as CIOS(1, A) instead of CIOS(A, 1), the formulas are slightly different but the result still holds.

Attacking consecutive CIOS steps

If Garner recombination is not used or the computation of CIOS(A, 1) is somehow protected against faults, a similar result can be achieved by forcing q ′ to zero in earlier calls to CIOS, provided that a certain number of successive CIOS executions are faulty.

Assuming that the values x and A in Montgomery representation are uniformly distributed modulo q before the first faulty CIOS, we show in the full version of this paper [START_REF] Fouque | Attacking RSA-CRT signatures with faults on Montgomery multiplication[END_REF]that faults across ℓ = ⌈log 2 ⌈log 2 q⌉⌉ iterations in the loop of the exponentiation algorithm are enough to ensure that S q will evaluate to zero with probability at least 1/2. For example, if q is a 512-bit prime, we have ℓ = 9. This means that forcing q ′ to zero in 9 iterations (from 9 to 18 calls to CIOS depending on the exponentiation algorithm under consideration and on the input bits) is enough to factor the modulus at least 50% of the time-and more faulty iterations translate to higher success rates. 1. Success rate of the null fault attack on consecutive CIOS steps, for a 512bit prime q and r = 16. 100 faulty signatures were computed for each parameter set.

For the Square-and-Multiply MSB and Montgomery Ladder algorithms, we compare success rates when faults start at the beginning of the loop vs. at a random iteration.

Simulation results. We have carried out a simulation of null faults on consecutive CIOS steps for each of the three exponentiation process algorithms, with varying numbers of faulty iterations; for the Square-and-Multiply MSB and the Montgomery Ladder algorithms, two sets of experiments have been conducted for each parameter set: one with faults starting from the first iteration, and another one with faults starting from a random iteration somewhere in the exponentiation loop. Results are collected in Table 1.

Constant Faults

In this section, we consider a different fault model, in which the fault attacker can force the variables u j in the CIOS algorithm to some (possibly unknown) constant value u.

Just as with null faults, we consider two scenarios: one in which the last CIOS computation is attacked, and another in which several inner consecutive CIOS computations in the exponentiation algorithm are targeted.

Attacking CIOS(A, 1)

Faults on all iterations. Consider first the case when faults are injected in all iterations of the very last CIOS computation. In other words, the device computes CIOS(A, 1), except that the variables u j , j = 0, . . . , k -1, are replaced by a fixed, possibly unknown value u. In that case, we show that a single faulty signature is enough to factor N and recover the secret key. The key result is as follows (the proof can be found in the full version [START_REF] Fouque | Attacking RSA-CRT signatures with faults on Montgomery multiplication[END_REF]).

Theorem 2. Let S be a faulty signature obtained in the fault model described above. Then, (2 r -1) • S is a close multiple of q with error size at most 2 r+1 , i.e. there exists an integer T such that:

(2 r -1) • ( S + 1) -qT ≤ 2 r+1 .
Thus, a single faulty signature yields a value V = (2 r -1) • ( S + 1) mod N which is very close to a multiple of q. It is easy to use this value to recover q itself. Several methods are available:

-If r is small (say 8 or 16), it may be easiest to just use exhaustive search: q is found among the values gcd(V + X, N ) for |X| ≤ 2 r+1 , and hence can be retrieved using around 2 r+2 GCD computations. -A more sophisticated option, which may be interesting for r = 32, is the baby step, giant step-like algorithm by Chen and Nguyen [START_REF] Chen | Faster algorithms for approximate common divisors[END_REF], which runs in time O(2 r/2 ). -Alternatively, for any r up to half of the size of q, one can use Howgrave-Graham's algorithm [START_REF] Howgrave-Graham | Approximate integer common divisors[END_REF] based on Coppersmith techniques. It is the fastest option unless r is very small (a simple implementation in Sageruns in about 1.5 ms on our standard desktop PC with a 512bit prime q for a any r up to ≈ 160 bits, whereas exhaustive search already takes over one second for r = 16).

Faults on most iterations. Howgrave-Graham's algorithm is especially relevant if the constant faults do not start at the very first iteration in the CIOS loop. More precisely, suppose that the fault attacker can force the variables u j to a constant value u not for all j but for j = j 0 , j 0 + 1, . . . , k -1 for some j 0 . Then, the same computation as in the proof of Theorem 2 yields the following bound on S q :

u • q 2 r -1 -2 rj 0 -2 < S q ≤ u • q 2 r -1 + 2 rj 0 + 1.
It follows that (2 r -1) • S is a close multiple of q with error size 2 r(j 0 +1) . Now note that Howgrave-Graham's algorithm [START_REF] Howgrave-Graham | Approximate integer common divisors[END_REF] will recover q given N and a close multiple with error size at most q 1/2-ε . This means that one faulty signature S is enough to factor N as long as j 0 + 1 < k/2, i.e. the constant faults start in the first half of the CIOS loop. Table 2. Success rate of the constant fault attack on successive CIOS steps, when using Square-and-Multiply LSB exponentiation with random 512-bit primes q and r = 16.

Attacking other CIOS steps

As in §3.2, if Garner recombination is not used or CIOS(A, 1) is protected against faults, we can adapt the previous attack to target earlier calls to CIOS and still reveal the factorization of N . However, the attack requires two faulty signatures with the same constant fault u. Details are given in the full version [START_REF] Fouque | Attacking RSA-CRT signatures with faults on Montgomery multiplication[END_REF].

In short, depending on the ratios q/2 ⌈log 2 q⌉ and u/(2 r -1), two faulty signatures S, S ′ with the same faulty value u have a certain probability of being equal modulo q. Thus, we recover q as gcd(N, S -S ′ ). This attack works with the Square-and-Multiply LSB and Montgomery Ladder algorithms, but not with Square-and-Multiply MSB exponentiation.

Simulation results are presented in Table 2. For various 512-bit primes q, the attack has been carried out for 1000 pairs of random messages, with a random constant fault u for each pair. It is successful if the two resulting faulty signatures S, S ′ satisfy gcd(N, S -S ′ ) = q.

Zero High-Order Bits Faults

In this section, we consider yet another fault model, in which the fault attacker targets the very last iteration in the evaluation of CIOS(A, 1) during the computation of S q . We assume that the attacker is able to force a certain number h of the highest-order bits of u k-1 to zero, possibly but not necessarily all of them (i.e. 1 ≤ h ≤ r). Then, while a single faulty signature is typically not sufficient to factor the modulus, multiple such signatures will be enough if h is not too small. More precisely, we prove the following theorem in the full version of this paper [START_REF] Fouque | Attacking RSA-CRT signatures with faults on Montgomery multiplication[END_REF]: Theorem 3. Let S be a faulty signature obtained in this fault model. Then, S is a close multiple of q with error size at most 2 -h • q + 1, i.e. there exists an integer T such that | S -qT | ≤ 2 -h • q + 1. Now, recovering q from faulty signatures of the form S is a partial approximate common divisor (PACD) problem, as we know one exact multiple of q, namely N , and several close multiples, namely the faulty signatures. Since the error size ≈ q/2 h is rather large relative to q, the state-of-the-art algorithm to recover q in that case is the one proposed by Cohn and Heninger [START_REF] Cohn | Approximate common divisors via lattices[END_REF] using multivariate Coppersmith techniques.

The algorithm by Cohn and Heninger is likely to recover the common divisor q ≈ N 1/2 given ℓ close multiples S (1) , . . . , S (ℓ) provided that the error size is significantly less than N (1/2) 1+1/ℓ . Hence, if the faults cancel the top h bits of u k-1 , we need ℓ of them to factor the modulus, where:

ℓ - 1 log 2 1 -h log 2 q . ( 3
)
In practice, if a few more faults can be collected, it is probably preferable to simply use the linear case of the Cohn-Heninger attack (the case t = k = 1 in their paper [START_REF] Cohn | Approximate common divisors via lattices[END_REF]), since it is much easier to implement (as it requires only linear algebra rather than Gröbner bases) and involves lattice reduction in a lattice of small dimension that is straightforward to construct. We examine this method in more details in the full version of this paper [START_REF] Fouque | Attacking RSA-CRT signatures with faults on Montgomery multiplication[END_REF], and find that it makes it possible to factor N provided that:

ℓ log 2 q h (4)
which is always a worse bound than (3) but usually not by a very large margin. Table 3 gives the theoretical number of faulty signatures required to factor N for various values of h, both in the general attack by Cohn and Heninger and in the simplified linear case. We carried out a simulation of the linear version of the attack on a 1024-bit modulus N with various values of h, and found that it works very well in practice with a number of faulty signatures consistent with the theoretical minimum. The results are collected in Table 4. The attack is also quite fast: a naive implementation in Sage runs in a fraction of a second on a standard PC.

Fault Models

In this section we discuss how realistic the setup of the attacks described above can be. In principle, all the RSA-CRT implementations using Montgomery multiplication may be vulnerable, but we have to note that the fault setup (and how realistic it is) depends heavily on implementation choices, since many variations around the algorithm from Figure 2 have been proposed in recent literature. After a discussion about the tools needed to get the desired effects, we focus on several implementation proposals [START_REF] Tenca | A scalable architecture for Montgomery multiplication[END_REF][START_REF] Mcivor | Modified Montgomery modular multiplication and RSA exponentiation techniques[END_REF][START_REF] Huang | An optimized hardware architecture for the Montgomery multiplication algorithm[END_REF][START_REF] Nozaki | Implementation of RSA algorithm based on RNS Montgomery multiplication[END_REF][START_REF] Suzuki | How to maximize the potential of FPGA resources for modular exponentiation[END_REF][START_REF] Mentens | Efficient pipelining for modular multiplication architectures in prime fields[END_REF][START_REF] Chow | A Karatsuba-based Montgomery multiplier[END_REF], chosen for their relevance.

Characteristics of the perturbation tool

First all the perturbations needed to carry out our attacks need to be controlled and local to some gates of the chip. Therefore, the attacker needs to identify the localization of the vulnerable gates and registers. The null fault attacks described in §3 need either a q ′ value set to 0, or multiple consecutive faults in line 6 of the main loop of CIOS(A, 1) or during multiple consecutive CIOS. The attacks described in §4 also need these multiple consecutive faults. Considering that state-of-art secure microcontrollers embed desynchronization countermeasures such as clock jitters and idle cycles, if the target of the perturbation is some shared logic with other treatments (like in the ALU of a CPU), the fault must be accurately space and time controlled, and the effects must be repeatable as well. Identification of the good cycles to inject the perturbation may be a very difficult task, and our attacks seem to be irrelevant. The only exception may be the null fault of §3, if the fault is injected when the q ′ register is loaded.

Nevertheless, many secure microcontrollers embed an isolated modular arithmetic acceleration coprocessor. A large proportion of them specifically use the Montgomery multiplication CIOS algorithm (or one of its described variants [START_REF] Koç | Analyzing and comparing Montgomery multiplication algorithms[END_REF]). Therefore, if the q ′ or the u j value is isolated in a specific small size register, a unique long duration perturbation can be sufficient for our attack to succeed. The duration of the perturbation varies with the implementation choices and can vary from one cycle to log 2 q, which does not exceed a hundred microseconds on actual chips. To get this kind of effect, laser diodes are the best-suited tool, since the duration of the spot is completely controlled by the attacker [START_REF] Skorobogatov | Optical fault masking attacks[END_REF].

Analysis of classical implementations of the Montgomery multiplication

The public Montgomery architectures can be divided in 3 different categories :

the first one [START_REF] Tenca | A scalable architecture for Montgomery multiplication[END_REF][START_REF] Mcivor | Modified Montgomery modular multiplication and RSA exponentiation techniques[END_REF][START_REF] Huang | An optimized hardware architecture for the Montgomery multiplication algorithm[END_REF] contains variations on the Tenca and Koç Multiple Word Radix-2 Montgomery Multiplication algorithm (MWR2MM) [START_REF] Tenca | A scalable architecture for Montgomery multiplication[END_REF], which can be seen as a CIOS algorithm with r = 1. The characteristic of these implementations is that they use no multiplier. They are then suited for constrained area. -the second category [START_REF] Suzuki | How to maximize the potential of FPGA resources for modular exponentiation[END_REF][START_REF] Mentens | Efficient pipelining for modular multiplication architectures in prime fields[END_REF] is an intermediate where r is a classical size for embedded architecture, such as 8,16 or 32 bits. They can be used for intermediate area/latency trade-offs. -the last category [START_REF] Nozaki | Implementation of RSA algorithm based on RNS Montgomery multiplication[END_REF][START_REF] Chow | A Karatsuba-based Montgomery multiplier[END_REF] propose a version of CIOS/SOS with only one loop, implying that r ≥ ⌈log 2 q⌉. The main difficulty of these implementation techniques is to deal with the very large multiplications they require . For that purpose they use interpolation techniques, like Karatsuba in [START_REF] Chow | A Karatsuba-based Montgomery multiplier[END_REF] or RNS in [START_REF] Nozaki | Implementation of RSA algorithm based on RNS Montgomery multiplication[END_REF]. These implementations are designed to achieve the shortest latency.

Architectures based on MWR2MM (r = 1). In this kind of architecture, q ′ cannot be manipulated, since it is always equal to 1, so no wire or register carries its value. On the other hand, the value of u j is computed at every loop of the CIOS, and since it is only one bit, a simple shot on the logic driving the register during the final multiplication CIOS(A, 1) is sufficient to get an exploitable result (u j = 0 corresponds to the null fault of §3, and u j = 1 to the constant fault of §4). The first proposal [START_REF] Tenca | A scalable architecture for Montgomery multiplication[END_REF] is a fully systolic5 array of processing elements (PE) executing consecutively line 6 of the CIOS algorithm in one cycle, and line 7 in k cycles from LSB to MSB. Figure 4 proposes an overview of the architecture. Each PE consists of a w-word carry save adder, able to compute a w word addition and to keep the carry for the next cycle. In the figure, T (j) stands for the j-th least significant w word of T .

At each clock cycle, the PE presents the computed result a i (j) to the next one, and the value u i is kept in the PE for the computation of the next word a i (j + 1). The value of u i is computed before the word a i (0) is presented, and then is kept in each PE during the whole computation of a i in a register. This architecture has the great advantage of being completely scalable (whatever the number of PEs and the size of M , this architecture can compute the expected result as long as the RAM are correctly dimensioned).

To achieve our attack, the register keeping u i can be the targeted, but every PE must be targeted simultaneously in order to get the correct result. Therefore it is more interesting to target the control logic responsible for the sequencing of the register loading, since all the PEs are connected.

In [START_REF] Mcivor | Modified Montgomery modular multiplication and RSA exponentiation techniques[END_REF], the authors manage to get rid of the CS to binary converter by redesigning the CS adder of every PE. The vulnerability to our attack is therefore the same, since the redesign does not affect the targeted area.

Huang et al. [START_REF] Huang | An optimized hardware architecture for the Montgomery multiplication algorithm[END_REF] proposed a new version of the data dependency in the MWR2MM algorithm and rearranged the architecture of [START_REF] Tenca | A scalable architecture for Montgomery multiplication[END_REF], in a semi systolic form. Figure 5 gives an overview of the architecture. In this architecture, the intermediate value a i is manipulated in carry save format ui-2 ui-j y(0) q(0) y( 1) q(1) y( 2) q(2) y(j) q(j)

Vulnerable areas

1 bit Right shift register containing x 1 bit Right shift register containing U • • • x(i) PE1 PEj PE2 PE0 x(i -1) x(i -2) x(i -j) register logic Combinational xi Vulnerable areas si(1) ci(1) ui y(0) q(0) ui si(1) ci(1) ci-1(2) si-1(2) • • •
A specific PE, PE 0 is specialized in generating the u i values at each cycle. while the j-th PE is in charge of computing the sequence a i (j). This architecture is very vulnerable to our attacks, since a simple ncycle long shot on the right logic in the PE 0 (see Figure 5) is sufficient to get the expected result.

According to the authors, the design works at 100 MHz on their target platform (a Xilinx Virtex II FPGA), therefore the duration of the perturbation is at least 10 µs for a 1024 bits multiplication (2048 bits RSA) if the Garner recombination is used (using the attack from §3.1 or §4.1). If classical CRT reconstruction is used, according to Table 1, 200 µs will be enough for a null fault.

As a conclusion we can see that this kind of implementation is very vulnerable, since the setup of the attack is quite simple.

High radix architecture (1 < r < ⌈log 2 q⌉). In this type of implementation the value q ′ = -q -1 mod 2 r is computed in a r-bit register, unless the quotient pipelining approach [23] is used.

For example, the implementation of [START_REF] Mentens | Efficient pipelining for modular multiplication architectures in prime fields[END_REF] is described in Figure 6. It relies on the coordinated usage of multiplier blocks of the Xilinx Virtex II together with specifically designed carry save adders. The values u j can be the target of any fault described in this paper, but it may be easier to put once for all the q ′ register to 0, with a 100% success rate for the attack if properly carried out. Another implementation is mentioned in [START_REF] Mentens | Efficient pipelining for modular multiplication architectures in prime fields[END_REF] with a four-deep pipeline, but it suffers from the same vulnerability.

The attack may be more difficult to achieve on the architecture of [28, Figure 4]. First, it uses quotient determination [START_REF] Orup | Simplifying quotient determination in high-radix modular multiplication[END_REF], and therefore does not need to store q ′ anywhere. Second, the multiplier in charge of computing x • y(i) u j is shared for all the Montgomery computation. In order to carry out the attack of §4 on this architecture, the attacker has to determine the specific cycles where u j is computed to generate a perturbation. For that particular design, the attacks seem out of reach.

Full radix architecture (r ≥ ⌈log 2 q⌉). In this kind of implementation, a single round is enough to compute the Montgomery algorithm. This implementation choice reports all the complexity on the design of a log 2 q × log 2 q multiplier. To reduce the full complexity of the big multiplication, interpolation techniques are used. In [START_REF] Chow | A Karatsuba-based Montgomery multiplier[END_REF], a classical nested Karatsuba multiplication is used, whereas [START_REF] Nozaki | Implementation of RSA algorithm based on RNS Montgomery multiplication[END_REF] proposes RNS. In these architectures, a specific laser shot must swap all the u 0 or q ′ at the same time to produce a null fault. To have a chance, a better solution is to use non invasive attacks (in the sense of [START_REF] Skorobogatov | Optical fault induction attacks[END_REF]), such as power or clock glitches. Indeed u 0 or q ′ are fully manipulated on the same clock cycle (or in very few), therefore it may be more practical to make the sequencer miss an instruction instead of aiming directly at the registers.

The zero high-order bits fault attack from §5 is also an option. In the architecture of [START_REF] Chow | A Karatsuba-based Montgomery multiplier[END_REF], the most significant bits of u 0 can be set to 0.

Conclusion

In this paper, we have shown that specific realistic faults can defeat unprotected RSA-CRT signatures with any padding scheme, probabilistic or not. While it is not difficult to devise suitable countermeasures (for example, checking that S q is not too small before outputting a signature is enough to thwart all of our attacks), this underscores the fact that relying on probabilistic signature schemes does not, in itself, protect against faults.

Fig. 2 .

 2 Fig.2. The Montgomery multiplication algorithm. The xi's and yi's are the digits of x and y in base b; q ′ = -q -1 mod b is precomputed. The returned value is (xy • b -k mod q). Since b = 2 r , the division is a bit shift.

Fig. 4 .

 4 Fig. 4. Systolic Montgomery Multiplier of [29] and potential target of the fault

Fig. 5 .

 5 Fig.5. Overview of the[START_REF] Huang | An optimized hardware architecture for the Montgomery multiplication algorithm[END_REF] architecture and potential target of the fault

Fig. 6 .

 6 Fig.6. Overview of the[START_REF] Mentens | Efficient pipelining for modular multiplication architectures in prime fields[END_REF] architecture and potential target of the fault

Table 3 .

 3 Theoretical minimum number ℓ of zero higher-order h-bit faulty signatures required to factor a balanced 1024-bit RSA modulus N using the general Cohn-Heninger attack or the simplified linear one.

	Number h of zero top bits			48 40 32 24 16
	Minimum ℓ with the general attack		8 9 11 15 22
	Minimum ℓ with the linear attack		11 13 16 22 32
	Number ℓ of faulty signatures 11 12 13 14 15 16 17 18
	Success rate with h = 48 (%) 23 100 100 100 100 100 100 100
	Success rate with h = 40 (%)	0	0	2 100 100 100 100 100
	Success rate with h = 32 (%)	0	0	0	0	0	0 99 100
	Average CPU time (ms)	33 35 38 41 45 49 54 59

Table 4 .

 4 Experimental success rate of the simplified (linear) Cohn-Heninger attack with ℓ faulty signatures when N is a balanced 1024-bit RSA modulus. Timings are given for our Sage implementation on a single core of a Core 2 CPU at 3 GHz.

Meaning that all the PEs are the same.