Mohammad Al Khatib

Antoine Girard
email: antoine.girard@l2s.centralesupelec.fr

Thao Dang
email: thao.dang@imag.fr

Mohammad Al Khatib
email: mohammad.alkhatib@l2s.centralesupelec.fr

Verification and Synthesis of Timing Contracts for Embedded Controllers *

Keywords: J.7 [Computer in other systems]: Command and control, Real time Verification, Design, Algorithms Stability, Reachability, Impulsive linear systems, Sampleddata systems

HAL is

INTRODUCTION

Physical systems equipped with embedded controllers are becoming pervasive (smart buildings, intelligent cars, drones, robots, etc.), thus increasing the need for high-confidence analysis and design tools that are able to handle tight interactions between the physical and digital worlds. In this context, contract-based approaches have been identified as a promising direction for cyber-physical systems design [START_REF] Sangiovanni-Vincentelli | Taming dr. frankenstein: Contract-based design for cyber-physical systems[END_REF]. For instance, for embedded controller implementation, [START_REF] Derler | Cyber-physical system design contracts[END_REF] proposed the use of timing contracts which specify the constraints on the time instants at which certain operations are performed such as sampling, actuation or computation. Under such contracts, the control engineers are responsible for designing a control law that is robust to all possible timing variation specified in the contract while the software engineers can focus on implementing the proposed control law so as to satisfy the timing contract. In this paper, we propose techniques that are useful within this framework. We first consider the problem of stability verification: given models of the physical plant and of the controller and a timing contract, verify that the resulting dynamical system is stable. We then tackle the problem of timing contract synthesis: given models of the physical plant and of the controller, determine a set of timing contracts that guarantee stability of the resulting system.

We adopt the impulsive linear dynamical system framework to model the overall system. Such systems form a class of hybrid systems which describe processes that evolve continuously and undergo instantaneous changes at discrete time instants. Applications of impulsive dynamical systems include sampled-data control systems [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], networked control systems [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF], multi-agent systems [START_REF] Bragagnolo | LMI sufficient conditions for the consensus of linear agents with nearly-periodic resets[END_REF], etc. In the present work, instantaneous changes occur at sampling and actuation times, t s k and t a k , k ∈ N, which are assumed to be non-deterministic. More precisely, we assume that some uncertainty lies in each of the k th sampling-to-actuation delay τ k = t a k -t s k and sampling period h k = t s k+1 -t s k . For the stability verification problem, we propose an approach based on the notion of reachable set. In the last decade, hybrid system reachability has had an important breakthrough in computing the reachable set corresponding to a linear continuous dynamics where the developed algorithms are based on representing the reachable sets by ellipsoids [START_REF] Kurzhanski | Ellipsoidal techniques for reachability analysis: internal approximation[END_REF][START_REF] Botchkarev | Verification of hybrid systems with linear differential inclusions using ellipsoidal approximations[END_REF], zonotopes [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF][START_REF] Althoff | Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes[END_REF] or by support functions [START_REF] Guernic | Reachability analysis of linear systems using support functions[END_REF][START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF]. In general, such algorithms handle time based switching by introducing auxiliary variables (clocks). In the following, we provide a specific approximation scheme for the reachable set at the sampling times to develop an effective stability verification approach. Then, we use this approach to tackle the timing contract synthesis problem. We propose a re-parametrization which provides some mono-tonicity property to the problem and allows us to develop an effective synthesis method based on guided sampling of the timing parameter space.

In all, we contribute in enriching our stability verification approach initiated in [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF] by stating the proofs of the necessary and sufficient theoretical conditions for the stability of the impulsive linear system and by considering more complex timing contracts which require dedicated reachability algorithms for stability verification and more involved techniques for timing contract synthesis.

The paper is organized as follows. First, some preliminary notations are defined before formulating the stability verification and timing contract synthesis problems in Section 2. Stability conditions are provided in Section 3. Section 4 presents the reachable set approximation scheme and an algorithm for stability verification. In Section 5, we propose a solution to the timing contract synthesis problem. In Section 6, examples, some of which are used to compare our results with existing ones, are then discussed before concluding our work.

Notations. Let R, R + 0 , R + , R - 0 , R -, N, N + denote the sets of reals, nonnegative reals, positive reals, nonpositive reals, negative reals, nonnegative integers and positive integers, respectively. For I ⊆ R + 0 , let NI = N∩I. Given a real matrix A ∈ R n×n , |A| is the matrix whose elements are the absolute values of the elements of A. Given S ⊆ R n and a real matrix A ∈ R n×n , the set AS = {x ∈ R n : (∃y ∈ S : x = Ay)}; for a ∈ R, aS = (aIn)S where In is the n × n identity matrix. The interior of S is denoted by int(S). The convex hull of S is denoted by ch(S). The interval hull of S is the smallest n-dimensional interval containing the set S and is denoted by 2(S). The symmetric interval hull of S is the smallest symmetric (with respect to 0) n-dimensional interval containing S and is denoted by (S). Given S, S ⊆ R n , the Minkowski sum of S and S is S ⊕ S = {x + x : x ∈ S, x ∈ S }. A polytope P is the intersection of a finite number of closed half-spaces, that is P = {x ∈ R n : Hx ≤ b} where H ∈ R m×n , b ∈ R m and the vector of inequalities is interpreted component-wise. Let Hi, i ∈ N [1,m] denote the row vectors of H, then if 0 ∈ int(ch({H1, . . . , Hm})), then P is compact. Given a template matrix H ∈ R m×n and a compact set S ⊆ R n , let us define the polytope ΓH (S) = {x ∈ R n : Hx ≤ b} where bi = maxx∈S Hix, i ∈ N [1,m] . In other words, ΓH (S) is the smallest polytope whose facets directions are given by H and containing S. We denote the set of all subsets of R n by 2 R n . We denote by K(R n ) the set of compact subsets of R n and by K0(R n ) the set of compact subsets of R n containing 0 in their interior. For x ∈ R, x is the smallest integer not less than x, and for c, c ∈ R n , c ≤ c if and only if ci ≤ c i , i = 1, . . . , n.

PROBLEM FORMULATION 2.1 Timing contracts for embedded control

In this work, we consider embedded control systems given under the form of general linear sampled-data control systems that take into account the sequences of sampling and actuation instants (t s k ) k∈N and (t a k ) k∈N :

ż(t) = Az(t) + Bu(t), ∀t ∈ R + u(t) = Kz(t s k ), t a k < t ≤ t a k+1 (1)
where z(t) ∈ R p is the state of the system, u(t) ∈ R m is the control input, and k ∈ N. For t ∈ [0, t a 0 ], u(t) can be any constant value in R m . We assume that the sequence of sampling and actuation instants (t s k ) and (t a k ) satisfy a timing contract given by

t s 0 = 0, τ k = t a k -t s k ∈ [τ , τ ], h k = t s k+1 -t s k ∈ [max(h, τ k ), h], k ∈ N (2)
where τ ∈ R + 0 , τ ∈ R + 0 , h ∈ R + , and h ∈ R + provide bounds on the sampling-to-actuation delays (which includes time for computation of the control law) and sampling periods provided that t s k ≤ t a k ≤ t s k+1 for all k ∈ N. Note that we impose h = 0 to prevent Zeno behavior. Moreover, these parameters must belong to the following set so that the time intervals given in (2) are always non-empty: 2) is a general timing contract which includes or over-approximates the different contracts introduced in [START_REF] Derler | Cyber-physical system design contracts[END_REF]. Their relation to the timing contract (2) is described as follows:

C = (τ , τ , h, h) ∈ R + 0 × R + 0 × R + × R + : τ ≤ τ ≤ h, h ≤ h Contract (
1. ZET Contract: The Zero Execution Time contract is given by ( 2) with τ = τ = 0 and h = h = h ∈ R + . In other words, the contract states that the sampling and actuation instants are periodic and simultaneous such that t s k = t a k = kh for k ∈ N. As mentioned in [START_REF] Derler | Cyber-physical system design contracts[END_REF], this contract is hardly achievable in practice since computation always takes time in between the sampling and actuation instants.

LET Contract:

The Logical Execution Time contract is given by ( 2)

with τ = τ = h = h = h ∈ R + .
The contract states that the sampling and actuation instants are periodic such that t s 0 = 0 and t s k = t a k-1 = kh for k ∈ N + .

DET Contract:

The Deadline Execution Time contract is given by (2) with τ = 0 and h = h = h ∈ R + . The contract states that the sampling instants are periodic, or t s k = kh for k ∈ N, and actuation instants are at some point t a k in the interval [t s k , t s k + τ ], with τ ≤ h.

TOL Contract:

The Timing Tolerance contract is defined by a nominal sampling period h ∈ R + , nominal sampling to actuation delay τ ∈ R + 0 , and two jitters J h , J δ ∈ R + 0 with J τ ≤ τ and

J h + J τ + τ ≤ h, such that t s k ∈ [kh, kh + J h ] and t a k ∈ [t s k + τ -J τ , t s k + τ + J τ ],
for k ∈ N (refer to Figure 1). We cannot exactly model this contract using (2). However we can over-approximate it using (2) with τ = τ -J τ , τ = τ +J τ , h = h-J h , and h = h+J h . Thus stability of system (1) under this latter contract guarantees also its stability under the TOL contract. LMI -scheduling protocols [START_REF] Hetel | Stabilization of linear impulsive systems through a nearly-periodic reset[END_REF] LMI τ = τ = 0 controller synthesis [START_REF] Hetel | Discrete and intersample analysis of systems with aperiodic sampling[END_REF] LMI

J h k ∈ [0, J h ] and J τ k ∈ [-J τ , J τ ].
τ = τ = 0 - [29] SOS τ = τ = 0 - [13] Invariant sets τ = τ = 0 - [1]
Reachability analysis τ = τ = 0 stochastic timing uncertainty [START_REF] Liu | Networked control systems in the presence of scheduling protocols and communication delays[END_REF] time-delay systems LMI h = 0 τ k > h k ; scheduling protocols [START_REF] Gao | Stabilization of networked control systems via dynamic outputfeedback controllers[END_REF] LMI h = h, τ = 0 controller synthesis; quantization [START_REF] Liu | Stability of linear systems with general sawtooth delay[END_REF] LMI

τ = τ = 0 - [15] LMI h = τ = τ = 0 - [4]
hybrid systems SOS -nonlinear dynamics; scheduling protocols [START_REF] Heemels | Networked control systems with communication constraints: Tradeoffs between transmission intervals, delays and performance[END_REF] LMI τ = 0, h = 0 scheduling protocols It is noteworthy that system (1) has deterministic dynamics under any of the ZET or LET contracts, unlike the case of the DET or TOL contracts where at least one of the sampling-to-actuation delays or sampling period is timevarying. In our previous work [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF] we considered the special case of nearly periodic linear impulsive systems (NPILS) which are modeled by ( 1) and ( 2) with τ = τ = 0. In other words, the sampling and actuation instants are simultaneous or t s k = t a k , and the duration in between two successive sampling instants is bounded in the interval [h, h]. Therefore, it is clear that we are dealing in this work with a more general timing contract which is more complex than the simple NPILS case.

Reformulation using impulsive systems

In our analysis it is more practical to transform equation ( 1) into an impulsive system with two types of resets each referring to a sampling or actuation instant. Such a reformulation is convenient to develop stability conditions based on reachability analysis. The system is thus given by:

ẋ(t) = Acx(t), t = t s k , t = t a k x(t s+ k ) = Asx(t s k ) x(t a+ k ) = Aax(t a k ) (3)
where x(t) ∈ R n is the state of the system with n = p + 2m, (t s k ) and (t a k ) are given by (2), x(t + ) = lim τ →0,τ >0

x(t + τ ), and

Ac =   A 0 B 0 0 0 0 0 0   , As =   Ip 0 0 K 0 0 0 0 Im   , Aa =   Ip 0 0 0 Im 0 0 Im 0   , x(t) =   z(t) Kz(θ s (t)) u(t)   , (4)
with θ s (t) = t s k for t ∈ (t s k , t s k+1 ]. In this paper, we consider stability in the following sense: Definition 1. The system (2-3) is globally uniformly exponentially stable (GUES) if there exist λ ∈ R + and C ∈ R + such that, for all sequences (t s k ) k∈N and (t a k ) k∈N verifying (2) the solutions of (3) verify

x(t) ≤ Ce -λt x(0) , ∀t ∈ R + .
We are now interested in verifying stability of embedded control systems in the form given by ( 1) under one of the general timing contracts defined previously. It is noteworthy that we can easily show that system (1) under the ZET and LET contracts is stable if and only if the eigenvalues of the matrix e hAc AaAs and Aae hAc As are inside the unit circle respectively. As for the DET or TOL contracts, we have that stability of system ( 1) is guaranteed by the stability of (2-3) with an adequate choice of the timing contract parameters. Consequently, in this work, we consider the following problem:

Problem 1 (Stability verification). Given Ac, As, Aa ∈ R n×n , (τ , τ , h, h) ∈ C, verify that (2-3) is GUES.
This problem will be considered in the following section. Afterwards, we shall consider the problem of synthesizing timing contract parameters that guarantee the stability of the system. Given bounds on the parameters 0 ≤ τmin ≤ τmax, 0 < hmin ≤ hmax, with τmin ≤ hmin, τmax ≤ hmax, let D = [τmin, τmax] 2 ×[hmin, hmax] 2 , the problem is formalized as follows:

Problem 2 (Timing contract synthesis). Given Ac, As, Aa ∈ R n×n and D, synthesize a set C * ⊆ C ∩ D such that for all (τ , τ , h, h) ∈ C * , (2-3) is GUES.

Related work.

Several approaches are developed in the literature to solve instances of Problem 1. A non-exhaustive list is given in Table 1. From the modeling perspective, the problem can be tackled using difference inclusions, timedelay systems or hybrid systems. On the computational side, most of the approaches are based on semi-definite programming using either Linear Matrix Inequalities (LMI) or Sum Of Squares (SOS) formulations. This makes a clear distinction with our approach which relies on reachability analysis. Let us remark that only a few approaches [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF][START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF][START_REF] Bauer | Stability analysis of networked control systems: A sum of squares approach[END_REF] appear to be able to address all instances of Problem 1. It is noticeable that [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF][START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF] have been implemented in the Networked Control Systems (NCS) toolbox [START_REF] Bauer | Networked control systems toolbox: Robust stability analysis made easy[END_REF] whose results will be compared to those of our approach. We should also acknowledge that some of these approaches are able to handle problems that we do not consider in the present work (possibility of having τ k > h k , controller synthesis, scheduling protocols, quantization, nonlinear dynamics, stochastic timing uncertainties). Finally, as far as we know, there is no available approach for addressing Problem 2 besides our preliminary work [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF] where we impose τ = τ = 0.

STABILITY CONDITIONS

In this section we state necessary and sufficient theoretical conditions for system (2-3) to be GUES. In addition, we derive practical sufficient conditions that can be used to develop an algorithm for solving Problem 1.

Necessary and sufficient conditions

Our stability conditions are based on the notion of reachable set defined as follows:

Definition 2. Given a continuous-time dynamical system ẋ(t) = Ax(t), t ∈ R + 0 , x(t) ∈ R n the reachable set on [t, t ] ⊆ R + 0 from the set S ⊆ R n is R A [t,t ] (S) = τ ∈[t,t ] e τ A S. (5)
We also define the map:

Φ : 2 R n → 2 R n , given for all S ⊆ R n by Φ(S) = τ ∈[τ ,τ ] w∈[max(0,h-τ ),h-τ ] e wAc Aae τ Ac AsS (6)
It is easy to see that if S is compact then so is Φ(S). It is clear that for two sets S, S ⊆ R n and a ∈ R, we have Φ(S ∪ S ) = Φ(S) ∪ Φ(S ) and Φ(aS) = aΦ(S).

The interpretation of Φ is as follows. If S is the set of all states that are reachable by (2-3) at time t s k then Φ(S) is the set of reachable states at time t s k+1 . We define the iterations of Φ as Φ 0 (S) = S for all S ⊆ R n , and Φ k+1 = Φ • Φ k for all k ∈ N. Then, for all k ∈ N, Φ k (S) is the set of reachable states by (2-3) at time t s k for initial states belonging to S. Next, we state the theoretical conditions on the stability of system (2-3) in terms of the map Φ. Similar results have been stated in [START_REF] Athanasopoulos | Alternative stability conditions for switched discrete time linear systems[END_REF] for discrete-time switched systems and in [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF] for NPILS.

Theorem 1. Let S ∈ K0(R n ), the following statements are equivalent:

(a) System (2-3) is GUES, (b) There exists a triplet (k, j, ρ) ∈ N + × N [0,k-1] × (0, 1) such that Φ k (S) ⊆ ρΦ j (S), (c) There exists a pair (k, ρ) ∈ N + ×(0, 1) such that Φ k (S) ⊆ ρ k-1 j=0 Φ j (S).
Proof. It is obvious that (b) =⇒ (c). Hence, it is sufficient to prove that (a) =⇒ (b) and (c) =⇒ (a).

(a) =⇒ (b): We will prove that there exists (k, ρ) ∈ N + × (0, 1) such that Φ k (S) ⊆ ρS. This implies (b) with j = 0. Let x(0) ∈ S, then Φ k (S) represents all the possible values of x(t s k ). Since (2-3) is GUES, there exist C ∈ R + and λ ∈ R + such that

x(t s k ) ≤ Ce -λt s k x(0) ≤ Ce -λkh x(0)
which can be rewritten as

Φ k (S) ⊆ Ce -λkh x(0) B ( 7
)
where B is the unit ball. Since S ∈ K0(R n ), then there exist c ∈ R + , c ∈ R + such that cB ⊆ S ⊆ cB. Then, [START_REF] Botchkarev | Verification of hybrid systems with linear differential inclusions using ellipsoidal approximations[END_REF] and

x(0) ∈ S give Φ k (S) ⊆ Ce -λkh cB ⊆ Ce -λkh c c S.
For k sufficiently large, Ce -λkh c < c and therefore (b) holds. (c) =⇒ (a): Let γ = ρ

1 k ; since ρ ∈ (0, 1) then for all j ∈ N [0,k-1] , ρ ≤ γ k-j and Φ k (S) ⊆ ρ k-1 j=0 Φ j (S) ⊆ k-1 j=0 γ k-j Φ j (S). ( 8
)
Let S = k-1 j=0 γ -j Φ j (S), then using properties of Φ:

Φ(S ) = Φ k-1 j=0 γ -j Φ j (S) = k-1 j=0 γ -j Φ j+1 (S) = k-2 j=0 γ -j Φ j+1 (S) ∪ γ -k+1 Φ k (S)
Making a change of index in the union and using (8) yield

Φ(S ) ⊆ k-1 j=1 γ -j+1 Φ j (S) ∪ γ -k+1 k-1 j=0 γ k-j Φ j (S) ⊆ γ k-1 j=0 γ -j Φ j (S) = γS . (9)
Since S ∈ K0(R n ), then S is compact. Moreover, since it contains S, then S ∈ K0(R n ). Then, there exist c ∈ R + , c ∈ R + such that c B ⊆ S ⊆ c B. Now consider a trajectory x of (2-3), then x(0) ∈ x(0) B ⊆ x(0) c S and ( 9) gives for all i ∈ N

x(t s i ) ∈ Φ i x(0) c S ⊆ x(0) γ i c S ⊆ x(0) γ i c c B.
In other words, it holds for all i ∈ N,

x(t s i ) ≤ γ i c c x(0) . Now, let t ∈ R + , let i ∈ N be such that t ∈ (t s i , t s i+1 ], then t -t s i ≤ h. Moreover, if t ∈ (t s i , t a i ], then x(t) ≤ e Ac h As γ i c c x(0)
and if t ∈ (t a i , t s i+1 ], then

x(t) ≤ e Ac h Aa As γ i c c x(0) .
In addition, we have i ≥ t/h and since γ ∈ (0, 1) it follows that for all t ∈ R +

x(t) ≤ e Ac h max( Aa , 1) As c c γ (t/h) x(0) ≤ e Ac h max( Aa , 1) As c c e ln(γ) h t x(0) .

Since γ ∈ (0, 1), (2-3) is GUES.

Sufficient conditions

The map Φ involved in Theorem 1 is in general impossible to compute exactly. Then, we may use an over-approximation Φ : K(R n ) → K(R n ) satisfying the following assumption:

Assumption 1. For all S ∈ K(R n ), Φ(S) ⊆ Φ(S).
We compute the map Φ instead of Φ in order to derive the practical condition on stability used in the stability verification algorithm later on to solve Problem 1. Section 4 discusses on the effective computation of the map Φ. We now derive sufficient conditions for stability based on Φ.

Corollary 1. Under Assumption 1, if there exist a set S ∈ K0(R n ) and a triplet (k, i, ρ) ∈ N + × N [0,k-1] × (0, 1) such that Φ k (S) ⊆ ρΦ i (S), then system (2-3) is GUES. Proof. Φ k (S) ⊆ ρΦ i (S) ⊆ ρ k-1 j=0 Φ j (S)
. Then similar to the second part of the proof of Theorem 1, let

S = k-1 j=0 γ -j Φ j (S) where γ = ρ 1 k . Then Φ(S ) = Φ k-1 j=0 γ -j Φ j (S) = k-1 j=0 γ -j Φ(Φ j (S)) ⊆ k-1 j=0 γ -j Φ(Φ j (S)) = k-1 j=0 γ -j Φ j+1 (S).
Then, following the same steps as in ( 9), we can show that Φ(S ) ⊆ γS . Following the same line as in the proof of Theorem 1, one concludes that (2-3) is GUES.

The previous corollary provides the background for designing a solution to Problem 1 in the next section.

OVER-APPROXIMATION SCHEME AND STABILITY VERIFICATION

In this section, we present an approach for computing an over-approximation of Φ. Furthermore, we develop an algorithm providing a solution to Problem 1.

Over-approximation

We first state the following result from [START_REF] Guernic | Reachability analysis of hybrid systems with linear continuous dynamics[END_REF] which gives an over-approximation scheme for the reachable set given by [START_REF] Bauer | Networked control systems toolbox: Robust stability analysis made easy[END_REF].

Theorem 2. [23] Let T ∈ R + , A ∈ R n×n , S ∈ K(R n ) and N ∈ N + , let R A [0,T ] (S) = N i=1 R A [(i-1)δ,iδ] (S)
where δ = T /N is the time step, and R A [(i-1)δ,iδ] (S) is defined by the recurrence equation:

R A [0,δ] (S) = ch(S, e δA S) ⊕ 1/4 δ (S), (10)
R A [iδ,(i+1)δ] (S) = e δA R A [(i-1)δ,iδ] (S), i ∈ N [1,N -1] with δ (S) = 2(|A| -1 (e δ|A| -I) (A(I -e δA )S)) ⊕ 2(|A| -2 (e δ|A| -I -δ|A|) (A 2 e δA S)). Then, R A [(i-1)δ,iδ] (S) ⊆ R A [(i-1)δ,iδ] (S), for all i ∈ N [1,N ] and R A [0,T ] (S) ⊆ R A [0,T ] (S).
The previous theorem can serve to compute an over-approximation of Φ. Indeed, from [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced lyapunov functions[END_REF], one can easily check that

Φ(S) ⊆ R Ac [max(0,h-τ ),h-τ ] AaR Ac [τ ,τ ] (AsS) ⊆ e max(0,h-τ )Ac R Ac [0,min(h-τ ,h-h+τ -τ )]
Aae τ Ac R Ac [0,τ -τ ] (AsS) [START_REF] Derler | Cyber-physical system design contracts[END_REF] with in turn can easily be over-approximated using the result of Theorem 2. In the case of NPILS, the previous inclusion becomes an equality. This is the approach followed in our previous work [START_REF] Khatib | Stability verification of nearly periodic impulsive linear systems using reachability analysis[END_REF]. However, for the general timing contract (2), the coupling in the timing uncertainties w and τ in ( 6) is totally disregarded in [START_REF] Derler | Cyber-physical system design contracts[END_REF] and therefore leads to conservatism. Therefore, in this paper, to reduce conservatism, we present a specific approximation scheme for Φ, that takes into consideration the coupling in the timing uncertainties.

It is based on the following result:

Lemma 1. Let S ∈ K(R n ), let N1, N2 ∈ N + , then Φ(S) ⊆ N 1 j 1 =1 n 2 (j 1 ) j 2 =1 e (θ(j 1 )+(j 2 -1)δ 2 )Ac R Ac [0,δ 2 ]
Aae (τ +(j 1 -1)δ 1 )Ac R Ac [0,δ 1 ] (AsS) [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF] where for j1

∈ N [1,N 1 ] δ1 = (τ -τ )/N1 δ2 = min(h -τ , h -h + δ1)/N2 θ(j1) = max(0, h -τ -j1δ1) n2(j1) = min(h -τ -(j1 -1)δ1, h -h + δ1)/δ2 . (13)
Proof. From (6), it follows that

Φ(S) = N 1 j 1 =1 τ ∈[τ +(j 1 -1)δ 1 ,τ +j 1 δ 1 ] w∈[max(0,h-τ ),h-τ ] e wAc Aae τ Ac AsS ⊆ N 1 j 1 =1 R Ac [θ(j 1 ),h-τ -(j 1 -1)δ 1 ] AaR Ac [τ +(j 1 -1)δ 1 ,τ +j 1 δ 1 ] (AsS) ⊆ N 1 j 1 =1 e θ(j 1 )Ac R Ac [0,h-τ -(j 1 -1)δ 1 -θ(j 1 )] Aae (τ +(j 1 -1)δ 1 )Ac R Ac [0,δ 1 ] (AsS) .
Remarking that

h -τ -(j1 -1)δ1 -θ(j1) = min(h -τ -(j1 -1)δ1, h -h + δ1) one gets Φ(S) ⊆ N 1 j 1 =1 n 2 (j 1 ) j 2 =1 e θ(j 1 )Ac R Ac [(j 2 -1)δ 2 ,j 2 δ 2 ]
Aae (τ +(j 1 -1)δ 1 )Ac R Ac [0,δ 1 ] (AsS) which leads to [START_REF] Donkers | Stability analysis of networked control systems using a switched linear systems approach[END_REF].

Remark 1. N1 and N2 are parameters used to discretize time intervals. For N1 = N2 = 1, the over-approximation given by ( 12) is the same as the one in [START_REF] Derler | Cyber-physical system design contracts[END_REF].

We now present our over-approximation scheme for Φ:

Theorem 3. Let S ∈ K(R n ), N1, N2 ∈ N + , and H ∈ R m×n , such that 0 ∈ int(ch({H1, . . . , Hm})), let Φ : K(R n ) → K(R n ) be given by Φ(S) = ΓH   N 1 j 1 =1 n 2 (j 1 ) j 2 =1
e (θ(j 1 )+(j 2 -1)δ 2 )Ac Φj 1 (S)

  where for j1 ∈ N [1,N 1 ] , Φj 1 (S) = R Ac [0,δ 2 ] Aae (τ +(j 1 -1)δ 1 )Ac R Ac [0,δ 1 ] (AsS)
with δ1, δ2, θ(j1), n2(j1) given by ( 13) , and R

Ac [0,δ 1 ] , R Ac [0,δ 2 ]
computed as in [START_REF] Cloosterman | Controller synthesis for networked control systems[END_REF]. Then, Φ(S) ⊆ Φ(S).

Proof. The proof is straightforward from Theorem 2 and Lemma 1.

Remark 2. In the previous result, the operation ΓH is not necessary to guarantee over-approximation of Φ. On the other hand, without this operation, the over-approximation of Φ would be given by the union of possibly numerous sets which may be quite impractical for subsequent manipulations. For that reason, this union is over-approximated by the smallest enclosing polytope whose facets direction are given by a matrix H. Moreover, if S is a polytope, then using the properties of support functions [START_REF] Guernic | Reachability analysis of linear systems using support functions[END_REF], the computation of Φ(S) reduces to solving a set of linear programs.

We illustrate the tightness of our new approximation scheme using system (16) (see Section 6) with the timing contract given by τ = 0, τ = 0.4, h = 0.2, h = 1.2. We consider a polytope S defined by a matrix H with 44 rows. Figure 2 shows sampled points (in grey) from Φ(S). The white polytope corresponds to the over-approximation Φ(S) given in Theorem 3 with N1 = 20 and N2 = 50. The black polytope is given by ( 11) over-approximated using Theorem 2 where R Ac [0,τ -τ ] and R Ac [0,min(h-τ ,h-h+τ -τ )] are computed with N = 20 and N = 50 respectively. One can check that the over-approximation given by Theorem 3 is quite tight and much less conservative than that given by (11).

Stability verification algorithm

We now present our stability verification algorithm which consists of two main steps: an initialization step where an initial set P0 is computed and a main loop which tries to verify the sufficient stability condition given in Corollary 1 by iterating the map Φ given by Theorem 3 from the set P0.

Initialization

The choice of the initial set P0 is crucial as it may impact significantly the number of iterations of Φ that are necessary to check the condition of Corollary 1. Intuitively, in order to minimize this number of iterations, P0 should be already close to an invariant set. Indeed, if Φ(P0) ⊆ P0, the stability condition holds after only one iterate of Φ. One way to choose P0 close to an invariant set is to define P0 as a common contracting polytope to L ∈ N + linear discrete-time systems, such that

∀j ∈ N [1,L] , e (h j -τ j )Ac Aae τ j Ac AsP0 ⊆ int(P0),
where the couples (τj, hj) satisfy timing contract (2) for all j ∈ N [1,L] . Then, P0 can be computed either using a backward iterative method as in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced lyapunov functions[END_REF] and [START_REF] Fiacchini | Set theory conditions for stability of linear impulsive systems[END_REF] or using a forward iterative method as in [START_REF] Athanasopoulos | Alternative stability conditions for switched discrete time linear systems[END_REF]. We denote the function computing P0 by init(Ac, A d , As, τ , τ , h, h, L). Then, P0 = {x ∈ R n : Hx ≤ b0}. The matrix H defining P0 is used in the main loop of the algorithm in the computation of the map Φ.

Main loop

The initial set is propagated using the map Φ given by Theorem 3. Then if the stability condition given by Corollary 1 is verified, system (2-3) is GUES and the algorithm returns true. Otherwise, if a maximum number of iterations, kmax, is reached then the algorithm fails to prove stability and returns unknown. The algorithm that solves Problem 1 is given as follows: Note that all polytopes P k are of the form P k = {x ∈ R n : Hx ≤ b k }, then the inclusion test at line 4 only consists in checking b k ≤ bi. Although Algorithm 1 is only based on sufficient conditions for the stability of system (2-3), its effectiveness will be demonstrated on numerical examples in Section 6.

Algorithm 1. Stability verification function is GUES(Ac,Aa,As,τ ,τ ,h,h) input: Ac, Aa, As ∈ R n×n , (τ , τ , h, h) ∈ C output: true if system (2-3) is proved GUES, unknown oth- erwise parameter: N1,N2,L, kmax ∈ N + 1: P0:=init(Ac, A d ,

TIMING CONTRACT SYNTHESIS

In this section, we propose a solution to Problem 2. We first define a re-parametrization of the timing-contract such that stability of system (2-3) becomes monotone with respect to the new parameters. Monotonicity is a very attractive property for designing efficient heuristics for timing contract synthesis since stability is preserved when the parameter values increase. This allows us to tackle the timing contract synthesis by sampling the parameter space.

Re-parametrization

Let us denote the vector of timing contract parameters α = (τ , τ , h, h) ∈ D = [τmin, τmax] 2 × [hmin, hmax] 2 , where 0 ≤ τmin ≤ τmax, 0 < hmin ≤ hmax, τmin ≤ hmin, τmax ≤ hmax. For α ∈ C ∩ D we denote the property:

Stab(α) ≡ (2-3) is GUES with parameters α.
Solving Problem 2 is equivalent to computing (a subset of) the set Co defined by

Co = {α ∈ C ∩ D : Stab(α)}.

Let us define a new parameter

β = (β1, β2, β3, β4) ∈ D where D = [τmin, τmax] × [-τmax, -τmin] × [hmin, hmax] × [-hmax, -hmin] and the map f : D → D such that f (β) = α = (τ , τ , h, h) where τ = β1, τ = min(-β2, -β4), h = β3, h = -β4.
We define the following constraint set for the parameter β:

C = β ∈ R + 0 × R - 0 × R + × R -: β1 ≤ min(-β2, -β4) β3 ≤ -β4 .
The following result holds:

Lemma 2. Let C o be given by

C o = {β ∈ C ∩ D : Stab(f (β))}. Then, f (C ∩ D ) = C ∩ D and f (C o ) = Co. Proof. Let us first show that f (C ∩ D ) ⊆ C ∩ D and f (C o ) ⊆ Co. Let β ∈ C ∩ D and α = f (β) = (τ , τ , h, h). Then, β ∈ D implies that α ∈ D, using the fact that τmin ≤ hmin, τmax ≤ hmax. Also, β ∈ C implies that τ ≤ τ and h ≤ h. Moreover, τ = min(-β2, -β4) ≤ -β4 = h. Hence, α ∈ C. Thus, α ∈ C ∩D. Moreover, if β ∈ C o then β ∈ C ∩D and Stab(f (β)) gives α ∈ C ∩ D and Stab(α). Thus, α ∈ Co. We now show that C ∩ D ⊆ f (C ∩ D ) and Co ⊆ f (C o ). Let α = (τ , τ , h, h) ∈ C ∩ D and let β = (τ , -τ , h, -h). Then, f (β) = (τ , min(τ , h), h, h). Since α ∈ C, it follows that min(τ , h) = τ and f (β) = α. Moreover, it is straightforward to verify that α ∈ C ∩ D implies β ∈ C ∩ D and that α ∈ Co implies β ∈ C o .
The previous result has two important implications. The first one is that the proposed re-parametrization does not introduce any conservatism in the solution to Problem 2 since the set Co of admissible parameters α can be obtained by computing the set C o of admissible parameters β, despite the fact that the map f is not injective nor surjective. The second one is stated in the following lemma:

Lemma 3. Let C * ⊆ C o , then C * = f (C * ) is a solution to Problem 2. Proof. It holds that C * = f (C * ) ⊆ f (C o ) = Co.
We further define the following set

E o = β ∈ D : (β / ∈ C ) ∨ ((β ∈ C ) ∧ Stab(f (β)))
One can easily check that the following relation holds:

C o = C ∩ E o . (14)
Hence, from the previous equality and Lemma 3, we can solve Problem 2 by computing (a subset of) the set E o . Moreover, E o satisfies the following monotonicity property:

Proposition 1. For all β, β ∈ D , the following implications hold:

(β ≤ β ) ∧ (β ∈ E o ) =⇒ β ∈ E o . (β ≤ β ) ∧ (β / ∈ E o ) =⇒ β / ∈ E o .
Proof. Let us assume β ≤ β and β ∈ E o . There are two cases:

1. If β / ∈ C , then either -β 4 ≤ -β4 < β3 ≤ β 3 , or -β 2 ≤ -β2 < β1 ≤ β 1 , or -β 4 ≤ -β4 < β1 ≤ β 1 . In all three cases β / ∈ C and therefore β ∈ E o . 2. If β ∈ C and Stab(f (β)), then either β / ∈ C which implies β ∈ E o , or β ∈ C . In this latter case, α = (τ , τ , h, h) = f (β) and α = (τ , τ , h , h ) = f (β ) sat- isfy α ∈ C, α ∈ C and τ ≥ τ , τ ≤ τ , h ≥ h, h ≤ h. (15)
It is straightforward to

check that if (2-3) is GUES for (τ , τ , h, h) ∈ C then (2-3) is GUES for all (τ , τ , h , h ) ∈ C satisfying (15). Thus, Stab(f (β )) holds and β ∈ E o .
This proves the first implication. For the second implication, it is sufficient to check that

(β ≤ β ) ∧ (β ∈ E o ) =⇒ β ∈ E o ≡ ¬(β ≤ β ) ∨ (β / ∈ E o ) ∨ (β ∈ E o ) ≡ (β ≤ β ) ∧ (β / ∈ E o ) =⇒ β / ∈ E o .
The previous property is instrumental for computing a subset of E o since it allows us to state the following theorem:

Theorem 4. Let β 1 , . . . , β M 1 ∈ E o , and β 1 , . . . , β M 2 ∈ D \ E o and let E = M 1 j=1 {β ∈ D : β j ≤ β}, E = D \ M 2 j=1 {β ∈ D : β ≤ β j }.
Then,

E ⊆ E o ⊆ E . Moreover, C * = f (C ∩ E ) is a solution to Problem 2 and Co ⊆ f (C ∩ E ).
Proof. E ⊆ E o ⊆ E is a direct consequence of Proposition 1. Then, from [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems[END_REF] and Lemmas 2 and 3, it follows that C * is a solution to Problem 2 and Co ⊆ f (C ∩ E ).

Timing contract synthesis algorithm

The previous theorem shows that it is possible to compute under and over-approximations of the set E o by sampling the parameter space D . In this section, we use this property to design a synthesis algorithm. Similar algorithms have been used in [START_REF] Legriel | Approximating the pareto front of multi-criteria optimization problems[END_REF][START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF] for computing an approximation of the Pareto front of a monotone multi-criteria optimization problem. Indeed, this latter problem can be tackled by computing an under and over-approximation of a set satisfying a monotonicity property similar to that of Proposition Lines 9 to 14 describe the main loop of the timing contract synthesis algorithm. At any time of the execution, E ⊆ E o ⊆ E holds. We pick a sample β ∈ E \ E which is the unexplored parameter region lying in the over-approximation of E o but not in its under-approximation. If β ∈ E o (or if β / ∈ E o ), then we update under-approximation E (or the over-approximation E ) according to Theorem 4. The algorithm stops when the Hausdorff distance between the E and E becomes smaller than ε. Of course, the choice of the sample β ∈ E \ E , at line 10, is crucial for the efficiency of the algorithm. In our implementation of the algorithm, we use the selection criteria proposed in [START_REF] Legriel | Approximating the pareto front of multi-criteria optimization problems[END_REF] which consists in choosing the sample that will produce the fastest decrease of the Hausdorff distance d(E , E ). In [START_REF] Tendulkar | Mapping and Scheduling on Multi-core Processors using SMT Solvers[END_REF] an alternative selection criteria based on multiscale grid exploration was proposed.

Finally, it is important to note that Algorithm 2 needs testing if the samples β ∈ E o which require checking the condition Stab(f (β)). In our implementation, this is done using Algorithm 1. If it returns true, then we can consider that Stab(f (β)) holds. If it returns unknown, we treat the sample as if Stab(f (β)) is false. As a consequence, in practice it may be the case that E is not an over-approximation of E o . However, it always holds that E ⊆ E o and therefore the set returned by Algorithm 2 is always a valid solution to Problem 2. Note that the property Stab(f (β)) need not be checked using Algorithm 1 but one can use any of the algorithms mentioned in Table 1.

ILLUSTRATIVE EXAMPLES

In this section, we first compare our approach for stability verification to that implemented within the NCS toolbox [START_REF] Bauer | Networked control systems toolbox: Robust stability analysis made easy[END_REF]. Then, we show an application of the timing contract synthesis algorithm. We implemented Algorithm 1 and Algorithm 2 in Matlab using the Multi-Parametric Toolbox [START_REF] Herceg | Multi-Parametric Toolbox 3.0[END_REF]. All reported experiments are realized on a desktop with i7 4790 processor of frequency 3.6 GHz and a 8 GB RAM.

Stability Verification

We consider two systems taken from [START_REF] Briat | Convex conditions for robust stability analysis and stabilization of linear aperiodic impulsive and sampled-data systems under dwell-time constraints[END_REF], given by (1) with the following matrices:

A = 0 1 0 -0.1 , B = 0 0.1 , K = -3.75 -11.5 . (
) A = 0 1 -2 0.1 , B = 0 1 , K = 1 0 . 16
We consider the stability verification problem for these two 2-dimensional systems. First, we write the systems into 4-dimensional impulsive systems (3). Then, we apply Algorithm 1 to check stability of the impulsive system under several timing contracts. We compare our results to those obtained using the NCS toolbox [START_REF] Bauer | Networked control systems toolbox: Robust stability analysis made easy[END_REF] in Table 2. For the DET timing contract (τ = 0, h = h = h), we fix parameter h and report the maximal value of τ for which stability has been verified. For the timing contract that corresponds to NPILS (τ = τ = 0), we fix h and report the maximal value of h for which stability has been verified. Finally, for the general timing contract given by (2), we fix parameters τ , τ , h and report the maximal value of h for which stability has been verified. Note that we conducted extra experiments labelled "Algorithm 1 (exp1)" to compare the results in terms of CPU time after fixing the same parameters as those used with the NCS toolbox.

The experiments conducted using the NCS toolbox are done in a particular manner since it uses three different approximation methods to embed the timing uncertainty (Jordan Normal Form (JNF), Cayley Hamilton, and Gridding and Norm Bounding (GNB)): we search for the maximum value of the free timing parameter that guarantees stability by running experiments using the three approximation methods. Then we report the computation time for the experiment in which we obtained this bound. In case the maximum bound could be obtained by more than one experiment, we report the CPU time corresponding to the fastest in terms of computation. Stability for system ( 16) is guaranteed using the GNB approximation for the DET, NPILS and general contracts, with 50, 35, and 50 gridpoints respectively. As for system [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF], stability is guaranteed using the JNF approximation for all three contracts. Parameter setups used by Algorithm 1, for the different experiments, are summarized by Table 3. Note that for the NPILS contract, the parameter N1 has no effect. It is clear, for the two systems at hand, that our method gives better results than the NCS toolbox in terms of CPU time and tightness.

Contract Synthesis

We now consider the timing contract synthesis problem for system [START_REF] Girard | Reachability of uncertain linear systems using zonotopes[END_REF]. We rewrite the system in the form of impulsive system (3). We search for a set C * ⊆ C ∩ D, where D = [0, 1.16] 2 × [0.21, 2.02] 2 , such that for all (τ , τ , h, h)∈ C * the system (2-3) is guaranteed to be GUES. We set the parameter ε = 0.07, and apply Algorithm 2. The output of the latter is a set C * ⊂ R 4 . Figure 3 shows a 3D section of C * by setting τ = 0.2, and a 2D section by setting τ = 0.2 and h = 1. Algorithm 2 used 3094 samples in the 4 dimensional parameter space with a total computation time of TCP U = 250 minutes. Parameters of Algorithm 1 used in Algorithm 2 are L = 4, kmax = 5 and the numbers of time steps used for the over-approximation of the reachable set are N1 = 20 and N2 = 50.

CONCLUSION

In this work, we proposed useful tools for contract-based design of embedded control systems under the form of algorithms for stability verification and timing contract synthesis. These algorithms can be used by control and soft-ware engineers to derive requirements that must be met by the real-time implementation of a control law. The effectiveness of our approach has been shown on examples. As future work, it would be interesting to handle the problem of controller synthesis given a timing contract, and to cosynthesize the controller and the timing contract parameters. Also more work is required for the stability verification problem as long as the solutions at hand gives only sufficient conditions.

Figure 1 :

 1 Figure 1: Time variables included in a TOL contract.J h k ∈ [0, J h ] and J τ k ∈ [-J τ , J τ ].

Figure 2 :

 2 Figure 2: Sampled points of Φ(S) (in grey), overapproximation Φ(S) given by Theorem 3 and overapproximation of (11) computed using Theorem 2 (in black).

Algorithm 2

 2 computes an under-approximation E and an over-approximation E of the set E o by sampling iteratively the parameter space D . Lines 1 to 8 correspond to the initialization of these approximations by testing the lower bound βmin and the upper bound βmax of the set D . If βmin ∈ E o , then by Theorem 4, f (C ∩ D ) = C ∩ D is a solution to Problem 2. Note that in that case, all timing-contract parameters in C ∩ D guarantee the stability of (2-3). If βmin / ∈ E o , then D \ {βmin} is an over-approximation of E o . Similarly, if βmax / ∈ E o , then by Theorem 4, E o = ∅. Note that in that case, no timingcontract parameters in C ∩ D can guarantee the stability of (2-3). If βmax ∈ E o , then {βmax} is an under-approximation of E o .

= 1 Figure 3 :

 13 Figure 3: (top) Timing contract synthesis for system (17) in the (0.2, τ , h, h) space where the visualized section of C * is in the domain region C defined above the planes h ≥ τ and h ≥ h. (bottom) The section of C * in the (τ , h) plane such that τ = 0.2 and h = 1.

Table 1 :

 1 Methods that can solve instances of Problem 1 with description of the modeling and computational approaches, list of restrictions and possible extensions.

	Models	Algorithm	Restrictions	Extensions
	[10] difference inclusion	LMI	-	τ k > h k ; controller synthesis
	[12]

1 .

 1 Algorithm 2. Timing contract synthesis function TC Synth(Ac,Aa,As,D) input: Ac, Aa, As ∈ R n×n , D = [τmin, τmax] 2 ×[hmin, hmax] 2 output: C * ⊆ C ∩ D such that for all (τ , τ , h, h) ∈ C

			, (2-3)
		is GUES.
	parameter: ε ∈ R +
	1: if βmin = (τmin, -τmax, hmin, -hmax) ∈ E o then
	2:	return C ∩ D;
	3: else E := D \ {βmin};
	4: end if
	5: if βmax = (τmax, -τmin, hmax, -hmin) / ∈ E o then
	6:	return ∅;
	7: else E := {βmax};
	8: end if
	9: while d(E , E ) > ε do	main loop
	10:	Pick β ∈ E \ E ;	select next sample
	11:	if β ∈ E o then E := E ∪ {β ∈ D : β ≤ β };
	12:	else E := E \ {β ∈ D : β ≤ β};
	13:	end if
	14: end while
	15: return f (C ∩ E );

*

Table 2 :

 2 Results of Algorithm 1 for systems[START_REF] Gao | Stabilization of networked control systems via dynamic outputfeedback controllers[END_REF] and (17) under several timing contracts. TCP U is the computation time in seconds.DET (τ = 0, h = h = h)

					NPILS (τ = τ = 0)		General contract (2)
		τ	h	TCP U	h	h	TCP U τ	τ	h	h	TCP U
	System (16) NCS toolbox (GNB) 0.63 1	3.42	10 -3 1.7291	3.30	0 0.4 0.2 1.13	9.17
	Algorithm 1(exp1)	0.63 1	0.18	10 -3 1.7291	0.20	0 0.4 0.2 1.13	4.49
	Algorithm 1(exp2)	0.67 1	1.16	10 -3 1.7294	0.20	0 0.4 0.2 1.23	9.95
	System (17) NCS toolbox (JNF) 0.78 1	2.07	0.4	0.45	1.91	0 0.1 0.4 0.44	3.62
	Algorithm 1(exp1)	0.78 1	0.41	0.4	0.45	0.21	0 0.1 0.4 1.71	1.13
	Algorithm 1(exp2)	1	1	2.97	0.4	1.88	1.22	0 0.1 0.4 1.71	5.15

Table 3 :

 3 Parameter setup for Algorithm 1 for systems (16) and (17) under several timing contracts. DET (τ = 0, h = h = h) NPILS (τ = τ = 0) General contract[START_REF] Althoff | Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes[END_REF]

		kmax N1 N2	L	kmax	N2	L	kmax N1 N2 L
	System (16)(exp1)	30	30	1	2	30	100	2	30	10 10 4
	System (16)(exp2)	30	100	1	2	30	1000	2	30	20 50 4
	System (17)(exp1)	30	15	1	2	30	1	2	30	10	1	4
	System (17)(exp2)	30	150	1	2	30	100	2	30	20 60 4

* This work was supported by the Agence Nationale de la Recherche (COMPACS project ANR-13-BS03-0004).


Céline Blondeau
email: celine.blondeau@aalto.fi

Kaisa Nyberg
email: kaisa.nyberg@aalto.fi

On Distinct Known Plaintext Attacks

Keywords: multidimensional linear attack, zero-correlation linear attack, key-difference-invariant-bias attack, known plaintext, distinct known plaintext, statistical model

Among statistical attacks, we often make a distinction between attacks in the linear context, for which the knowledge of plaintexts and corresponding ciphertexts is enough to perform the attack, and attacks in the differential context, for which the plaintexts are chosen. Such attacks are usually referred as known-or chosen-plaintext attacks. It is commonly believed that attacks in the known-plaintext model are more practical than attacks in the chosen-plaintext model. Nevertheless, it is usual in the literature, to only compare the data, time and memory complexity of these attacks without considering the model. In this paper, we reconsider some known-plaintext attacks, by considering them in the distinct-known-plaintext model. We explain and develop the statistical model for the multiple zero-correlation linear cryptanalysis, multidimensional linear cryptanalysis, as well as for the key-differenceinvariant-bias related-key attack introduced at ASIACRYPT 2013. Based on these models validated by experiments, we improve attacks on some ciphers.

Introduction

Among statistical attacks, we often make a distinction between attacks in the linear context, for which the knowledge of plaintexts and corresponding ciphertexts is enough to perform the attack, and attacks in the differential context, for which the plaintexts are chosen. Such attacks are usually referred as known-or chosen-plaintext (KP, CP) attacks. It is commonly believed that attacks in the known-plaintext model are more practical than attacks in the chosen-plaintext model. Nevertheless, it is usual in the literature to only compare the data, time and memory complexity of these attacks without considering the model.

In [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF][START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF][START_REF] Bogdanov | Zero Correlation Linear Cryptanalysis with Reduced Data Complexity[END_REF] zero-correlation linear attacks were introduced. For these attacks, depending on the number of used approximations and on the relation between the involved linear masks, different statistical models are presented. For instance in [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF] the two models to compute respectively the data complexity of multiple zero-correlation linear attacks and multidimensional zero-correlation linear attacks are recalled. While the statistical model for multidimensional zerocorrelation linear attack assumes that the plaintexts involved in the attacks are distinct, the one for multiple zero-correlation linear attack assumes a normal distribution of the expected capacity for the wrong keys.

In [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF] key-difference-invariant-bias attacks are introduced. In these relatedkey attacks, the attacker takes advantage of linear approximations that have the same bias for keys with a specific difference. The statistical model [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF] for this related-key known-plaintext attack is similar to the one presented in [START_REF] Bogdanov | Zero Correlation Linear Cryptanalysis with Reduced Data Complexity[END_REF] in the context of multiple zero-correlation linear attack.

In this paper we develop on distinct-known-plaintext (DKP) attacks. In particular, we show the importance of avoiding repetition when the data complexity of the attack is close to the full codebook. More importantly, while it is commonly believed that there is a difference in the evaluation of the data complexity of a multiple zero-correlation linear attack and of a multidimensional zero-correlation linear attack, we show that the difference is only in the way the plaintexts are handled. Based on theoretical and experimental observations, we provide new formulas to compute the data complexity of multiple zero-correlation linear attacks and key-difference-invariant bias attacks. As an illustration, we have instantiated our formulas on existing attacks reducing their data complexity. A model for distinct-known-plaintext multiple/multidimensional linear attacks is also introduced.

The outline of this paper is as follows. The notation is introduced in Sect. 2. In Sect. 3 we present the statistical model which is validated by experiments in Sect. 4. In Sect. 5 we apply our model on some existing attacks and improve their complexity. Sect. 6 concludes this paper.

Preliminaries

Linear Attacks

While the idea of using distinct-known plaintexts can be extended to any statistical attack, we focus on the most common known-plaintext statistical attacks which are generalizations of linear cryptanalysis [START_REF] Matsui | Linear Cryptanalysis Method for DES Cipher[END_REF].

Given an n-bit permutation F , we denote by (u, v) ∈ F n 2 × F n 2 , a pair of input and output masks. In linear attacks, we take advantage of linear approximations of the form u • x ⊕ v • F (x) = 0. The strength of a linear relation is measured by its correlation. The correlation of a Boolean function

f u,v (x) = u • x ⊕ v • F (x) is defined as cor(u, v) = 2 -n # {x ∈ F n 2 |f u,v (x) = 0} -# {x ∈ F n 2 |f u,v (x) = 1} .
In [START_REF] Biryukov | On Multiple Linear Approximations[END_REF] the idea of taking advantage of multiple independent linear approximations is introduced. In [START_REF] Hermelin | Multidimensional Extension of Matsui's Algorithm 2[END_REF] multidimensional linear attacks are presented. In these more recent attacks, the attacker takes advantage of all linear approximations with linear masks (u, v) u = 0 in a linear space.

The capacity C is used to collect the information of a set of linear approximations. Given sets of input and output linear masks U and V , it is defined as the sum of the squared correlations:

C = u∈U,v∈V,u =0 cor 2 (u, v).
While multiple/multidimensional linear attacks take advantage of a set of linear approximations with large capacity, multiple and multidimensional zerocorrelation linear attacks [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF][START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF][START_REF] Bogdanov | Zero Correlation Linear Cryptanalysis with Reduced Data Complexity[END_REF] take advantage of linear approximations with correlation equal to zero. These attacks have been proven efficient on wordoriented structures such as Feistel-type ciphers. When multiple approximations with zero-correlation are used, the capacity C of the set of linear approximations is equal to zero.

In the remainder of this paper, we denote by the number of linear approximations involved in our attacks. Given s the dimension of the linear space U × V , in (zero-correlation) multidimensional linear attacks we have = 2 s -1. The block cipher size is denoted by n.

Statistics

The data complexity N of a statistical attack corresponds to the number of plaintexts necessary to perform the attack. In general, we want to find the encryption key (right key) by differentiating the score of the right key from the one of the wrong keys. In (zero-correlation) multiple/multidimensional linear attacks the scoring function corresponds to the estimated capacity of the multiple/multidimensional linear approximations:

T = N • 1≤i≤ ( ĉ or i ) 2 ,
where ĉ or i is the empirical correlation of the i-th linear approximation. In (zerocorrelation) multidimensional linear attacks the computation of this score T can be simplified and is equivalent to:

T = j=0 (V [j] -N/( + 1)) 2 N/( + 1) ,
where V [j] corresponds to the number of occurrences of the j-th element of the multidimensional distribution. Following the notation of [START_REF] Aydin | On Probability of Success in Linear and Differential Cryptanalysis[END_REF], we denote by P s the success probability and by a the advantage of the attack where 2 -a is the proportion of discarded keys.

Throughout this paper, we denote by Φ the cumulative distribution function of the central normal distribution. To simplify the notation, we also introduce:

ϕ a = Φ -1 (1 -2 -a
) and ϕ P S = Φ -1 (P S ). Given µ R and σ R (resp. µ R and σ w ), the mean and variance of the normal random variable T R for the right key (resp. T W for the wrong keys), we have (see i.e. [START_REF] Aydin | On Probability of Success in Linear and Differential Cryptanalysis[END_REF]):

P S ≈ Φ |µ R -µ W | -σ W ϕ a σ R . ( 1
)
3 Statistical Models

Multiple and Multidimensional Zero-Correlation Linear Attacks

In [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF] we have the following two estimates of the data complexity of a multiple and multidimensional zero-correlation linear attacks derived from [START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF][START_REF] Bogdanov | Zero Correlation Linear Cryptanalysis with Reduced Data Complexity[END_REF].

Lemma 1. [7]
The number N of known plaintexts required in a multiple zerocorrelation linear cryptanalysis is:

N ≈ 2 n (ϕ P S + ϕ a ) /2 -ϕ a . (2)
The proof [START_REF] Bogdanov | Zero Correlation Linear Cryptanalysis with Reduced Data Complexity[END_REF] follows from (1) using the fact that the distribution of T R (resp T W ) can be estimated by a normal distribution with parameters µ R = N and

σ R = √ 2 N (resp. µ W = N + 2 n and σ W = √ 2 N + √ 2 2 n ). Lemma 2. [4, 6] 1
The number N of distinct-known plaintexts required in a multidimensional zero-correlation linear cryptanalysis is:

N ≈ 2 n (ϕ P S + ϕ a ) /2 + ϕ P S . (3)
The proof [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF][START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF] follows from the use of the hypergeometric distribution as given in a more general case in Th. 1. Assuming as in the proof of Lemma 1 that the correlation of the involved linear approximations are independent, we can adapt this result to the context of multiple zero-correlation linear cryptanalysis. In practice we observed, see Sect. 4.1, that the data complexity of a multiple zero-correlation linear attack can be estimated by (3) when assuming distinct-known plaintexts.

Corollary 1. Given the number of used linear approximations. The data complexity of a known-plaintext multiple/multidimensional zero-correlation linear attack is given by (2). If we consider distinct-known plaintexts, the data complexity is given by (3).

Since for most attacks 0.5 ≤ P S ≤ 0.99, meaning that 0 ≤ ϕ P S ≤ 2.4, the difference between (3) and ( 2) is particularly noticeable when /2 and ϕ a are in the same order of magnitude. From ( 3) and ( 2) we deduce that the success probability of a known-plaintext zero-correlation linear attack is:

P S ≈ Φ N 2 n /2 -ϕ a • ( N 2 n + 1) , (4)
and the one of a distinct-known-plaintext zero-correlation linear attack is:

P S ≈ Φ N /2 2 n -N -ϕ a 2 n 2 n -N . ( 5
)

Multidimensional Linear Attacks

In the previous section we show that the data complexity of zero-correlation linear attacks is reduced once we consider distinct-known plaintexts. In this section we focus on the classical multidimensional linear attack.

Lemma 3. [START_REF] Hermelin | Multidimensional Extension of Matsui's Algorithm 2[END_REF] In the known-plaintext model, the random variable T R involved in a multiple/multidimensional linear attack follows a normal distribution with parameters:

µ R ≈ + N • C, and (6)
σ 2 R ≈ 2( + 2 • N • C).
The random variable T W follows a normal distribution with parameters µ W = ,

σ 2 W = 2 .
From these results, the data complexity of a known-plaintext multidimensional linear attack is computed as [START_REF] Hermelin | Multidimensional Extension of Matsui's Algorithm 2[END_REF]:

N ≈ √ 4a + 4Φ -1 (2P S -1) 2 C . (7)
In the following we study the distribution of the variable T R in the context of a distinct-known-plaintext attack.

Theorem 1. Assuming distinct-known plaintexts the random variable T R involved in a multiple/multidimensional linear attack follows a normal distribution with parameters:

µ R ≈ (1 - N 2 n ) + N • C, and (8)
σ 2 R ≈ 2 (1 -N/2 n ) 2 + 4(1 -N/2 n )N • C.
Proof. For the purpose of this proof, we denote by E(X) the mean of a random variable X and by V ar(X) its variance. We denote by Z i the random variable corresponding to the number of solutions of the i-th equation of the form u • x ⊕ v • F (x) = 0 and by N • p i the expected number of solutions of this equation.

Assuming that the plaintexts are distinct, from the hypergeometric distribution,

we have E(Z i ) = N • p i = N 2 (1 + cor i ) and V ar(Z i ) = N p i (1 -p i ) = N 4 ( 1
-
cor 2 i ) 2 n -N 2 n -1 since p i = 1/2(1 + cor i ). Given X i = 2 • Z i /N -1 we deduce that E(X i ) = cor i and V ar(X i ) = 4 N 2 V ar(Z i ) ≈ 1 N • 2 n -N 2 n -1 . We have T R = N i X 2 i . We denote V = i≤ X 2 i V ar(X i ) . The random vari- able V follows a non-central χ 2 distribution with parameters E(V ) = + λ and V ar(V ) = 2( + 2λ) where λ = i E(X i ) 2 V ar(X i ) = N • C 2 n -1 2 n -N . From V = T R 2 n -1 2 n -N , we deduce that E(T R ) ≈ (1 -N/2 n ) + N • C and V ar(T R ) ≈ 2 (1 -N/2 n ) 2 + 4(1 -N/2 n ) • N • C.
Experiments which confirm the value µ R are presented in Sect. 4.2. From this result we can extract using C = 0 the distribution of the random variable T R (see [START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF]) in the case of a distinct-known-plaintext zero-correlation linear attack.

Experimental Results

As usually zero-correlation linear attacks and multidimensional linear attacks are not targeting the same cipher construction, we have implemented experiments on a Feistel-type cipher and on a SPN-type cipher. The used ciphers are depicted in Fig. 1 and could correspond to scaled versions of CLEFIA [START_REF] Shirai | The 128-Bit Block cipher CLEFIA (Extended Abstract)[END_REF] (a 16-bit type-II GFN with 4 branches) and PRESENT [START_REF] Bogdanov | PRESENT: An Ultra-Lightweight Block Cipher[END_REF] (SMALLPRESENT-[32], a 32-bit SPN-cipher).

F F s E e s E e F F s E e s E e F F s E e s E e 0 0 u 0 0 u 0 0 K 1 2 K 11 1 K 12 2 X 1 X 2 X 11 X 12 X 13
ZC on 9 rounds

⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ S 7 S 6 S 5 S 4 S 3 S 2 S 1 S 0
Fig. 1: Left: Description of the key-recovery attack done on a Type-II GFN. Right: One round of SMALLPRESENT- [START_REF] Boura | Scrutinizing and Improving Impossible Differential Attacks: Applications to CLEFIA, Camellia, LBlock and Simon[END_REF].

While in [START_REF] Soleimany | Zero-correlation linear cryptanalysis of reduced-round LBlock[END_REF] experiments showing the distribution of µ R and µ W have been presented, there is, to the best of our knowledge, no mentioning of experimental zero-correlation linear attacks, in the literature.

Zero-Correlation Attacks on a Type-II GFN with 4 Branches

The results of our experiments averaged over 1000 keys are provided in Fig. 2. In these graphics we compare the success probability of multidimensional and multiple zero-correlation linear attacks with the theoretical ones given by ( 5) for distinct plaintexts and by (4) for non-distinct plaintexts. This experiments confirm the theory given in Sect. 3.1 showing that the same formula can be used to compute the complexity of multiple zero-correlation and multidimensional zero-correlation linear attacks.

Fig. 2: Attacks on a type-II-GFN cipher. Left: multidimensional zero-correlation linear attacks, right: multiple zero-correlation linear attacks.

Experiments on SMALLPRESENT-[8]

In Fig. 3 we compare the experimental and theoretical mean µ R of the variable T R in the cases of distinct-known-plaintext and of known-plaintext distinguishing attacks. For this cipher, we observe that the theoretical value of µ R given in Th.1 is accurate. Exp. non-distinct Fig. 3: The mean µR of TR for a 8-bit multidimensional distribution ( = 2 8 -1) over 9 rounds of SMALLPRESENT- [START_REF] Boura | Scrutinizing and Improving Impossible Differential Attacks: Applications to CLEFIA, Camellia, LBlock and Simon[END_REF] with capacity C = 2 -21.29 (left), over 12 rounds with capacity C = 2 -24.13 ≤ /2 n (right).

Remark 1. Distinct-known-plaintext attacks should require less data than knownplaintext attack when the data complexity of the attack is close to the full codebook. In practice, as illustrated by [START_REF] Boura | Scrutinizing and Improving Impossible Differential Attacks: Applications to CLEFIA, Camellia, LBlock and Simon[END_REF], this occurs when the capacity of the multiple/multidimensional linear approximations is such that C ≤ /2 n . In the distinct-known-plaintext model this would correspond to the case where µ R decreases as the data complexity increases (see i.e. the right side of Fig. 3). An attack in this model will only be possible if we have a good estimate of the capacity. However, in practice, see for instance [START_REF] Joo | Linear Cryptanalysis of Reduced-Round PRESENT[END_REF], we only have an underestimate of the capacity which can be problematic when it comes to estimate the data complexity of a distinct-known-plaintext multidimensional linear attack.

Applications

Multiple Zero-Correlation Linear Attacks

As explained in this paper, by considering distinct-known plaintexts we can use [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF] to compute the data complexity of a multiple zero-correlation linear attack. As the data complexity of multidimensional linear attacks has already been computed under this setting, and because other comparable (in number of attacked rounds) attacks have been performed in the chosen-plaintext model, this should give us a better comparison factor. The result of our computation and a comparison with the best attacks on the block cipher Camellia [START_REF] Aoki | Camellia: A 128-Bit Block Cipher Suitable for Multiple Platforms -Design and Analysis[END_REF] are provided in Table 1. The attack is from [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF]. The data complexity has been computed using (3) instead of using (2) with the parameters of the attack chosen as P S = 0.85 and a = 96 or a = 160. The time complexity has been computed according to the description given in [START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF].

Similarly we can improve the data complexity of the multiple zero-correlation linear attack on CAST-128 [START_REF] Wen | General Application of FFT in Cryptanalysis and Improved Attack on CAST-256[END_REF]. The parameters of the attack being n = 128, = 64770, a = 50 and P S = 0.85, the data complexity of the attack using known plaintexts2 is N = 2 123.73 and the data complexity of the attack using distinct-known plaintexts is N = 2 The memory is expressed in number of bytes. #R denotes the number of attacked rounds.

Key-Difference-Invariant-Bias Attacks

To the best of our knowledge, the only paper presenting attacks in this context is the similar paper [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF]. In Table 2 we resume the complexity of the best related keyattacks on LBlock [START_REF] Wu | LBlock: A Lightweight Block Cipher[END_REF] and show that by assuming distinct-known plaintexts the data and time complexity of the attack can be improved. Similar improvement can be obtained for the related-key attack on TWINE presented in [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF]. This paper* Table 2: Best related-key attacks on LBlock. *: Computation of the time complexity according to the description given in Sect. 5.3 of [START_REF] Bogdanov | Key Difference Invariant Bias in Block Ciphers[END_REF].

Conclusion

In this paper, we reconsider the statistical model for multiple zero-correlation linear and key-difference-invariant-bias attacks. We show that when using distinct plaintexts, the attacks can be performed using less plaintexts. We also consider for the first time a statistical model for distinct-known-plaintext multiple/multidimensional linear attacks. Nevertheless questions remain regarding how this model will be useful when we only have an underestimate of the capacity.

Table 1 :

 1 123.67 . Correlation 2 14 160 85% 2 125.46 DKP 2 125.8 2 112 This paper Best key-recovery attacks on Camellia (attacks starting from the first round).

	Version	#R	Attack	a PS N	Time Mem.	Ref.
	Camellia-128 11	Impossible	--	1 2 118.4 CP	2 118.43 2 96.4	[8]
	Camellia-128 11 Zero-Correlation 2 14 96 85% 2 125.3 KP	2 125.8 2 112	[4]
	Camellia-128 11 Zero-Correlation 2 14 96 85% 2 125.1 DKP 2 125.8 2 112 This paper
	Camellia-192 12	Impossible	--	1 2 119.7 CP	2 161.06 2 147.7	[8]
	Camellia-192 12 Zero-Correlation 2 14 160 85% 2 125.7 KP	2 125.8 2 112	[4]
	Camellia-192 12 Zero-

  RKCP 2 78.3 2 61.4 [16] 24 Key Inv Bias 32 2 7.81 4.5 85% 2 62.29 RKKP 2 74.59 2 61 [15] 24 Key Inv Bias 32 2 7.81 8.5 85% 2 62.95 RKKP 2 70.67 2 61 [15] 24 Key Inv Bias 32 2 7.81 8.5 85% 2 62.38 RKDKP 2 70.67 2 61 This paper 24 Key Inv Bias 32 2 7.81 16 85% 2 62.84 RKDKP 2 66.57 2 61

	#R	Type	#Keys		a PS N	Time Mem.	Ref.
	23 Imp. Diff	4	-	100% 2 61.4

The distribution of the random variables has been derived in[START_REF] Bogdanov | Integral and Multidimensional Linear Distinguishers with Correlation Zero[END_REF], the estimation of the data complexity appears in[START_REF] Bogdanov | Zero-Correlation Linear Cryptanalysis with FFT and Improved Attacks on ISO Standards Camellia and CLEFIA[END_REF].

With these parameters, the data complexity can not be equal to 2 123.2 as given in[START_REF] Wen | General Application of FFT in Cryptanalysis and Improved Attack on CAST-256[END_REF].


Sondre Rønjom
email: sondrer@gmail.com

Powers of Subfield Polynomials, Cyclic Codes and Algebraic Attacks with Applications to the WG Stream Ciphers

In this paper we investigate univariate algebraic attacks on filter generators over extension fields Fq = F2n with focus on the Welch-Gong (WG) family of stream ciphers. Our main contribution is to break WG-5,WG-7,WG-8 and WG-16 by combining results on the so-called spectral immunity (minimum distance of certain cyclic codes) with properties of the WG type stream cipher construction.The spectral immunity is the univariate analog of algebraic immunity and instead of measuring degree of multiples of a multivariate polynomial, it measures the minimum number of nonzero coe cients of a multiple of a univariate polynomial. Based on the structure of the general WG-construction, we deduce better bounds for the spectral immunity and the univariate analog of algebraic attacks.

Introduction

There exist at least five published variants of the WG construction; WG-5 [START_REF] Aagaard | Hardware implementations of the WG-5 cipher for passive RFID tags[END_REF],WG-7 [START_REF] Luo | A Lightweight Stream Cipher WG-7 for RFID Encryption and Authentication[END_REF],WG-8 [START_REF] Fan | WG-8: A Lightweight Stream Cipher for Resource-Constrained Smart Devices, In Quality, Reliability, Security and Robustness in Heterogeneous Networks[END_REF],WG-16 [18] and WG-29 [START_REF] Nawaz | The WG Stream Cipher[END_REF]. In this section we present basic results that make up the machinery of these constructions, needed for our cryptanalysis of WG-ciphers in Section 3.

M-sequences, Unitary Sequences and Linear Complexity

For a much better introduction to the relationship between finite fields and sequences, the reader is referred to [START_REF] Lidl | Finite Fields In Encyclopedia of Mathematics and its Applications[END_REF] and [START_REF] Golomb | Signal Design for Good Correlation: For Wireless Communication, Cryptography, and Radar[END_REF]. Let F q n denote a n-th order extension of the binary field F q of q = 2 k elements, defined by a polynomial m(x) over F q . In order to simplify the presentation we assume that m(x) is a primitive polynomial. The polynomial m(x) defines a linear feedback shift register (LFSR) over F q of length n that generates a maximal sequence (or msequence) of period q n 1; if initzialised in a non-zero state the LFSR spans the coe cient vectors of exactly the elements of the multiplicative group F ⇤ q n . Moreover, if the LFSR is initizialised in a state S 0 = (s 0 , s 1 , . . . , s n 1 ) 2 F n q , the LFSR generates a sequence obeying the recurrence relation

s t+n =s t c 0 + s t+1 c 1 + . . . + s t+n 1 c n 1 . (1)
defined by the coe cients of m(x). The minimal polynomial of a periodic sequence s over F q n is the polynomial of least degree generating that sequence. The degree of this polynomial is what is called the linear complexity of the sequence. Let ↵ 2 F q n be, for sake of simplicity, a primitive element. For 2 F ⇤ q n and d 2 {0, 1, . . . , q n 1} we call the sequence b d,t = (↵ t ) d , t = 0, 1, 2, . . . over F q n a unitary sequence. Unitary sequences b d,t are the simplest forms of nonzero sequences in the sense that their linear complexity is 1, since their minimal polynomials are the linear polynomials x + ↵ d . It is well-known that the minimal polynomial of the sum of two sequences a t and b t is equal to the least common multiple of their individual minimal polynomials. Thus the sum of two unitary sequences a t = 1 (x↵ t ) d1 and b t = 2 (x↵ t ) d2 has simply minimal polynomial m(x) = (x + ↵ d1 )(x + ↵ d2 ). In general, if I is a random distinct subset of {0, 1, 2, . . . , q n 1} and c i are random nonzero constants of F q n , the polynomial

P (x) = X i2I c i x i defines a sum of unitary sequences z t = P i2I c i (x↵ t ) i with minimal polynomial m(x) = Q i2I (x + ↵ i ) and linear complexity |I|.
2 Filter Generators and Algebraic Attacks over F 2 n

Filter generators have been well-studied in literature and consists usually of a binary m-sequence generating LFSR of length n, a Boolean function f in k variables and a subset of tapping positions I ⇢ {i 1 , i 2 , . . . , i k } ⇢ {0, 1, 2, . . . , n}. In this section we quickly recapture the current state of algebraic attacks on such constructions, but in terms of univariate polynomial equations. In the rest of the paper, all operations on polynomials over F q n are modulo ), a binary keystream sequence can be generated by

x q n + x. Let L(x) = P n 1 i=0 x 2 i denote the trace from F q = F 2 n to F 2 and ↵ 2 F 2 n a
z t =f (s t+i1 , s t+i2 , . . . , s t+i k ) =f (L t+i1 (x), L t+i2 (x), . . . , L t+i k (x))
The bits of the sequence z t are successively exored with the plaintext bits to form a ciphertext sequence. The choice of LFSR, tapping positions and Boolean function all have various e↵ects on the cryptographic quality of the resulting keystream z t .

Algebraic Attacks

In algebraic cryptanalysis (see for instance [START_REF] Armknecht | Introducing a new variant of fast algebraic attacks and minimizing their successive data complexity[END_REF] and [START_REF] Hawkes | Rewriting variables: The complexity of fast algebraic attacks on stream ciphers[END_REF]) of a filter generator the adversary tries to solve an associated equation system relating unknown statevariables with keystream values. Then if the adversary has observed a sequence of keystream bits beginning at time t, (z t , z t+1 , . . . , z t+m ), she can set up a system of equations of the form nonzero coe cients (exactly those x i where wt(i)  d). This means that if the adversary observes D keystream bits, she can set up a system of at most D equations in D unknowns over F 2 n and solve using linear algebra. In multivariate cryptanalysis the complexity is given by O(D log2 (7) ). Notice that F t+i (x) = F t (x↵ i ) and that the coe cients of F t (x) = P wt(i)d c i x i ↵ ti span cyclic vectors of the form v t = (↵ ti1 , ↵ ti2 , . . . , ↵ tiD ). If we let D such vectors for t = 0, 1, 2, . . . , D 1 span a D ⇥ D matrix M , the resulting matrix is a Vandermonde matrix and can be manipulated more e ciently than generic matrices (the inverse can be computed in O(Dlog(D) 2 )). Moreover, if X = (c i1 x i1 , c i2 x i2 , . . . , c iD x iD ) then M • X = (z 0 , z 1 , . . . , z D 1 ) and M 1 (z 0 , z 1 , . . . , z D 1 ) = (x i0 , x i1 , . . . , x iD ) If we compute X from z, we can easily recover the initial state x from one of the equations c ij x ij = x ij . In practice one can pre-compute one of the columns of the inverse to recover x from a pre-chosen value x ij . This is essentially the improved algebraic attack presented in [START_REF] Rønjom | A New Attack on the Filter Generator[END_REF].

z t =f (L t+i1 (x), L t+i2 (x), . . . , L t+i k (x)) =F t (x), t = 0,

Algebraic Attacks and Low-Degree Polynomials

An often more keystream e cient method is to make use of low-degree multivariate multiples of f and f + 1. Moreover, if there exist a multivariate Boolean polynomial g in the ideal spanned by f over F 2 of lower degree e, the adversary can use the relation

g(S t )(z t + f (S t )) =0 which yields a new valid equation each time z t = 0 since the zeros of f is a subset of any multiple g. If we let G t (x) = g(L t+i1 (x), L t+i2 (x), . . . , L t+i k (x)), we can construct a system of equations G ti (x) = 0 for all t i when z ti = 0. Let T = {t 1 , t 2 , . . . , t E } where E = P e i=0 n i .
The equations involve at most E nonzero coe cients so we can set up a E ⇥ E matrix M spanned by coe cient vectors v tj = (↵ tj i1 , ↵ tj i2 , . . . , ↵ tj iE ) for t j 2 T and an unknown initial state related vector X = (c

i1 x i1 , c i2 x i2 , . . . , c iE x iE ) such that v ti • X = G t (x) = 0 for all t i 2 T .
The rank of the equation system in an "annihilator"-attack has been assumed to have almost full rank E in litererature, but it has been an open question. We can now resolve this question by noting that the matrix M is a generalized Vandermonde matrix and it was shown by Shparlinski [START_REF] Shparlinski | On the Singularity of Generalised Vandermonde Matrices over Finite Fields[END_REF] that almost all such matrices have full rank. The algebraic immunity of a Boolean function was introduced in [START_REF] Meier | Algebraic attacks and decomposition of Boolean functions[END_REF] and measures the resistance of a function against algebraic attacks. Moreover, the algebraic immunity, abbreviated AI(f ), is defined as the minimal degree of a multiple of either f or f + 1. It has been shown that AI(f ) for a k-variable function satisfy the bound 0  AI(f )  dk/2e. The adversary can therefore always reduce datacomplexity from P d i=0 n i to roughly 2

P d k
2 e i=0 n i if the degree of the function f is larger than AI(f ). But all hope is not lost even if the design employs a Boolean function with optimal algebraic immunity. It was shown in [START_REF] Courtois | Fast algebraic attacks on stream ciphers with linear feedback[END_REF], that if there exist polynomials g and h with deg(g) < deg(h) < deg(f ) where h = g • f , the adversary can instead set up an equation system of the form

h t (S 0 ) + g t (S 0 ) • z t = 0 for t = 0, 1, 2, . . .. Let e = deg(g), d = deg(h), E = P e i=0 n i
and

D = P d i=0 n i . Further, let H t (x) = h t (L t+i1 (x), L t+i2 (x), . . . , L t+i k (x)) and G t (x) = g t (L t+i1 (x), L t+i2 (x), . . . , L t+i k (x)) such that H t (x) + G t (x) • z t = 0.
The authors of [START_REF] Courtois | Fast algebraic attacks on stream ciphers with linear feedback, Advances in Cryptology-Crypto[END_REF] showed that if the adversary pre-computes the minimal polynomial m h (x) of the sequence b t = h(S t ) she can simply apply the recursion defined by m h (x) =

P D i=0 c i x i to the equation system D X i=0 c i (H t+i (x) + G t+i (x)z t+i ) = 0 for t = 0, 1, 2, . . . , E 1. The polynomial m h (x) is simply Q ci6 =0 (x + ↵ i ) where c i are the coe cients of H 0 (x)
where we assume that all the coe cients for terms x i of weight less or equal to d are nonzero. Since the sequence h(S t ) = H t (x) obeys the recursion defined by m h (x), the new equations become

D X i=0 c i (G t+i (S 0 )z t+i ) = 0
for t = 0, 1, 2, . . . which is now a system of equations involving the E coe cients of G t (x). The best total complexity for solving such equation systems has been shown to be O(EDlog 2 (D) + E log2 (7) ). It is assumed that one needs D + E keystream bits to solve this system, since the relation D is used to determine E equations. However, in practice one can compute a polynomial of degree D E and zeros ↵ i with e < wt(i)  d that will cancel only the terms of x i where i has weight larger than e, so only D keystream bits are needed in practice. However, O(D + E) = O(D) for typical applications, so it usually makes little or no di↵erence.

3 Filter Generators in the Spirit of Welch-Gong

The Welch-Gong type filter generator consists of a primitive LFSR over an extension field F q = F 2 k of length n and a Boolean function f (x) over F q . Let ↵ 2 F q n denote a root of the LFSR generator polynomial. The LFSR defines a q-ary sequence

s t = L(x↵ t ) where L(x) = T r q n /q (x) = P n 1 i=0 x q i denotes the trace from F q n to F q and x 2 F
q n is a random nonzero initial state. The WG-design applies a Boolean function f (x) to exactly one q-ary element L(x↵ t ) of the LFSR register, in e↵ect generating a binary keystream

z t = f (L(x↵ t
)) for t = 0, 1, 2 . . .. In the following section our focus will be on minimizing the complexity of univariate algebraic attacks on this particular construction.

Powers of Subfield Polynomials and Minimum Distance

When solving univariate equations we do not care so much about degree as we care about the number of nonzero coe cients in the polynomials. The equations we are interested in are of the form

z t =f (L(x↵ t )) = q 1 X i=0 c i L(x↵ t ) i =F (x↵ t )
where f is over F q and F (x) is over F q n . In the rest of the paper we will write capital letters F, G, H to represent functions f, g, h over F q composed with L(x), where L(x) will be fixed in the context. Notice that we need not take the composition f (L(x)) modulo x q n + x since the highest degree term possible in L(x) q 1 is q (n 1) (q 1) = q n q n 1 . To any polynomial f (x) over F q , define a weight enumerator polynomial

T f (x) = k X i=0 w i x i (2)
where w i counts the number of nonzero terms x d in f with exponent d of hamming weight i. We have that T f (1) is the usual hamming weight if the coe cients of f are binary. If f (x), L(x) are as above it is easy to determine the number of nonzero coe cients of their composition F (x).

Theorem 1. The number of nonzero coe cients of

F (x) = f (L(x)) is given by T f (n).
Proof. In the expansion of L(x) e with e = 2 u1 + 2 u2 + . . . , 2 u d with u i 1 < u i we get terms of the form

x q i 1 2 u 1 +q i 2 2 u 2 +•••q i d 2 u d . (3)
Each varying i and j in the range 0  i < n and 0  j < k, we get that q i 2 j corresponds to nk distinct powers 2 0 , 2 1 , . . . , Due to this special structure of F (x) we can improve the bounds on the so-called spectral immunity of univariate polynomials. Spectral immunity of a general Boolean polynomial F (x) over F q n was essentially defined in [START_REF] Gong | Fast linear subspace attacks on stream ciphers[END_REF] in terms of sequences, but here it is more convenient to use the definition provided by Helleseth et. al. [START_REF] Helleseth | Simplifying algebraic attacks with univariate analysis[END_REF] in terms of cyclic codes.

2 nk 1 . Moreover, since 0  u 1 < u 2 < . . . < u d < k,
choice of u = (u 1 , u 2 , . . . , u d ) (corresponding to di↵erent powers of L(x)). Let A u = {q i 2 u | 0  i  n 1}. Notice that A u \ A u 0 = ; for u 6 = u 0 . Each distinct choice of u gives rise to n d exponents A u1 + A u2 + . . . + A u d = {b 1 + b 2 + . . . + b d | b i 2 A ui }
Theorem 2. The spectral immunity of a Boolean function F (x) over F q n , denoted SI(F ), is equal to the minimum weight of a q-ary cyclic code generated by

G F (x) = gcd(F (x) + 1, x q n + x) or G F +1 (x) = gcd(F (x), x q n + x).
The spectral immunity is the univariate analog of algebraic immunity as it measures the least number of unknowns one needs to solve for in an algebraic attack. In a general algebraic attack over F 2 n (when the LFSR is defined over F

2 ), we have polynomials of the form

P (x) = f (L t+i1 (x), L t+i2 (x), . . . , L t+i k (x))
Since the algebraic immunity of f (x) as a multivariate polynomial in k variables is at most dk/2e, it follows that the spectral immunity of P (x) is upper-bounded by P dk/2e i=0 n i . Although univariate and multivariate attacks have similar complexity in general, as we have seen, the WG-type construction produces polynomials of a very special type that allows us to improve this bound significantly.

Lemma 1. Let F (x) be of the above form (essentially defining a WG-cipher). The minimum distance of the cyclic codes generated by G F and G F +1 over F q n is upper-bounded by

SI(F )  dk/2e X i=0 ✓ k i ◆ n i ( ✓ 2 • dk/2e 1 dk/2e ◆ 1)n dk/2e .
Proof. The proof is straight-forward. Assume the worst case, which is when f (x) is balanced. The matrix M i containing the

P dk/2e i=0 k i coe cient vectors of x d (mod g i (x
)) where g i (x) = gcd(f (x) + i, x q + x), i 2 F 2 and with d of hamming weight less or equal to dk/2e, has rank at most 2 k 1 . Consequently, the kernel K i of M i has dimension at least

P dk/2e i=0 k i 2 k 1 = 2•dk/2e 1 dk/2e
. Since It is not clear whether the upper-bound for T g (n) for a multiple g of f or f + 1 (the least number of coe cients of a multiple G of F or F + 1) is tight or not, and leave this as an open problem.

Minimizing Data Complexity In Attacks on WG Ciphers

In our attacks on the WG-ciphers we seek to minimize the data-complexity.

To do this we focus on relations f • g = h where h and g only have terms x i of weight at most dk/2e. Let M i denote the matrix spanned by the coe cient vectors of x d (mod r i (x)) where r i (x) = gcd(f (x) + i, x q + x) for i 2 F 2 and where wt(d)  dk/2e. When k is odd we have that 2•dk/2e 1 dk/2e = k dk/2e , and when k is even it is equal to k 1 dk/2e = k dk/2e /2. When k is odd the kernel K i has at least dimension k dk/2e , so there must exist multiples g 0 (x) of f (x) and g 1 (x) of f (x) + 1 with at most one (and the same) term x i with exponent weight dk/2e such that their sum g(x) = g 0 (x) + g 1 (x) has only terms of weight < dk/2e. In particular, since f (g 0 + g 1 ) ⌘ g 0 (mod x q + x) and (f + 1)(g 0 + g 1 ) ⌘ g 1 (mod x q + x) any equation

z t + f (L t (x)) = 0
with f over F 2 k with k odd can always be replaced by an equation 0

= g(L t (x))(z t + f (L t (x))) = g(L t (x))z t + g 0 (L t (x)) (4)
where G(x) = g(L(x)) has at most P dk/2e 1 i=0 k i n i nonzero coe cients while G 0 (x) has at most P dk/2e 1 i=0 k i n i + n dk/2e nonzero coe cients. Similar argument can be made when k is even. Since in this case the kernels have at least dimension k 1 dk/2e , we can find a polynomial g 0 in K 0 and g 1 in K

1 such that the two polynomials share at least k 1 dk/2e

1 terms and coe cients of weight dk/2e. Since their sum g(x) = g 0 (x) + g 1 (x) has at most k 1 dk/2e + 1 terms of weight dk/2e, we can construct equations of the form

0 = g(L t (x))(z t + f (L t (x))) = g(L t (x))z t + g 0 (L t (x)
) but where g 0 (L t (x)) now has at most

P dk/2e 1 i=0 k i n i +( k 1 dk/2e 1)n dk/2e nonzero coe cients. Let E = P dk/2e 1 i=0 k i n i and D = P dk/2e i=0 k i n i ( 2•dk/2e 1 dk/2e
1)n dk/2e such that the complexity of a fast algebraic attack on a generic WG cipher is upper-bounded by C = EDlog(D) + E log2 (7) . Then if k is small in comparison to n, C is roughly equal to O((

✓ k dk/2e 1 ◆ n dk/2e 1 ) log2(7) ) (5)
and data complexity of O(n dk/2e ) keystream bits.

Cryptanalaysis of the WG Family

The class of WG-ciphers are pure filter generator constructions consisting of an LFSR of length n over F q = F 2 k and a k-variabe Boolean function f (x) over F q . In this section we analyse WG-5, WG-7, WG-8 and WG-16. In the following let L(x) denote the trace polynomial from F q n to F q where n and q = 2 k is clear from the context.

Breaking WG-5

For WG-5 [START_REF] Aagaard | Hardware implementations of the WG-5 cipher for passive RFID tags[END_REF] we have n = 32 and F q = F 2 5 . The Boolean function is given by f (x) = T r(x d ) over F q where the specification leaves a choice of either using d = 7 or d = 15. Since the functions have optimal algebraic immunity, the designers state that the best possible algebraic attack has complexity 2 54 using 2 19 keystream bits. From the rough bound 5, there exist g and h where g(L(x)) has n 2 = 2 10 nonzero coe cients and h(L(x)) has n 3 = 2 15 coe cients. One can verify that g 0 (x) = x 24 + x 8 + x 7 + x 5 + y and g 1 (x) = x 24 + x 9 + x 8 + x 7 + x 5 + y satisfy f (g 0 + g 1 ) + g 0 = 0. Moreover, from f (x) • g(x) = g 0 (x), we get equations of the form

z t • L(x↵ t ) 9 = g 0 (L(x↵ t ))
for t = 0, 1, 2, . . . where the left-hand side involves n 2 = 2 10 unknowns and the right-hand side roughly n 3 = 2 15 unknowns. Thus the complexity of a fast algebraic attack is roughly (n 2 ) log(7,2) ⇡ 2 30 using n 3 = 2 15 keystream bits. We find 5 3 = 10 such relations that can be used to mount the same attack for both T r(x 7 ) and T r(x 15 ). Although using the function T r(x 15 ) result in a higher linear complexity than T r(x 7 ) against an algebraic attack, they behave the same against our attack.

4.2 Breaking WG-7 WG-7 [START_REF] Luo | A Lightweight Stream Cipher WG-7 for RFID Encryption and Authentication[END_REF] consists of an LFSR of length 23 over a field F 2 7 and a Boolean function f (x) = T r(x 3 + x 9 + x 21 + x 57 + x 87 ) corresponding to a multivariate Boolean function in seven variables. The linear complexity of the keystream generated by WG-7 is approximately 2 25.5 and the authors assume that an attacker has access to no more than 2 24 keystream bits.

If the function has optimal algebraic immunity the complexity of an algebraic attack is roughly 2 69 using 2 25 keystream bits. Using our rough bounds, the complexity is at most 2 38 using 2 18 keystream bits, which is much less than 2 25 . But the authors of [START_REF] Orumiehchiha | Cryptanalysis of WG-7: a lightweight stream cipher[END_REF] find that the algebraic immunity f (x) is in fact 3, and mount algebraic a fast algebraic attack in 2 28 using 2 19.38 keystream bits. Due to the low algebraic immunity, it is easy to find exactly one low weight multiple g 0 (x) of f (x) and g 1 (x) for f (x) + 1 where all coe cients for both polynomials are 1 and all exponents have hamming weight less or equal to 3. The sum of g 0 (x) and g 1 (x) is simply g(x) = g 0 (x) + g 1 (x) = T r(x) = P 7 i=0 x 2 i which can be used to construct a set of equations z t T r(L(x↵ t )) = g 0 (L(x↵ t )) for t = 0, 1, 2, . . . . The number of unknowns in the right-hand side equation is given by T g0 (n) = 2 13.84 ⇡ 2 14 and the left-hand side has n • k = 161 unknowns. Moreover, a FAA on this uses only 2 14 keystream bits and has complexity about

(n • k) ⇥ 2 14 log 2 (2 14 ) + (n • k) log(7) ⇡ 2 25
. This is a factor 2 3 faster than the algebraic attack of [START_REF] Orumiehchiha | Cryptanalysis of WG-7: a lightweight stream cipher[END_REF], but more importantly uses only a factor 2 5 of the their keystream complexity which has the practical significance in this setting.

Breaking WG-8

WG-8 [START_REF] Fan | WG-8: A Lightweight Stream Cipher for Resource-Constrained Smart Devices, In Quality, Reliability, Security and Robustness in Heterogeneous Networks[END_REF] consists of an LFSR of length 23 over F q = F 2 8 and apply the Boolean function f (x) = T r(x 9 + x 37 + x 63 + x 127 ) over F q . The authors claim that the best algebraic attack on this construction is in 2 69 using 2 26 keystream bits. To find good relations for this specification, we computed the kernel of M 0 and M 1 consisting of rows spanned by the coe cients of x t (mod g i (x)) for all t of weight  4 and where g i (x) = gcd(f (x) + i, x q + x).

Each of the kernels had minimal dimensions 35, and their sum K = K 0 + K 1 dimension 70. We simply row reduced the basis matrix for K to eliminate the high weight terms first and collected the last row of the matrix. It contained the coe cient vector of a polynomial s(x) = s 0 (x) + s 1 (x) with only one term x d of weight 4 and had E = T s (n) = 2 17 . Moreover, the polynomials s i (L(x)) have D = T si (n) ⇡ 2 22 nonzero coe cients. In an attack we can therefore use a relation z t s(L(x↵ t )) + s 0 (L(x)) = 0 for t = 0, 1, 2, . . . and obtain an attack with complexity EDlog(D)+E log(7) ⇡ 2 48 using 2 22 keystream bits. We found many relations that gives the same attack complexity.

4.4 Breaking WG-16 WG-16[18] consists of a primitive LFSR of length 32 over F q = F 2 16 and is meant for use in 4G. The function f (x) has multivariate degree 8 and optimal algebraic immunity. The authors claim that the best algebraic attack on this construction is in 2 159 using 2 58 keystream bits. Since the function is in an even number of variables, the minimal T g (n) for a multiple g of f and f + 1 satsify T g (n)  P dk/2e 1 i=0 k i n i +( k 1 dk/2e + 1)n dk/2e ⇡ 2 53 . A direct univariate algebraic attack has then complexity T g (n) log(7,2) = 2 148 using 2 54 keystream bits which is already less than the claimed bounds. Using the bounds for relations f • g = h, there must exist a polynomials g and h where g(x) has exactly one term x i with wt(i) >= dk/2e 2 and h(x) has terms of weight up to dk/2e+2. Thus we can set E = T g (n)  P dk/2e 2 i=0 k i n i + n dk/2e ⇡ 2 37 and D = T h (n)  P dk/2e+2 i=0 k i n i ⇡ 2 63 that yields an attack with computational complexity ED log(D)+E log(7,2) ⇡ 2 106 using 2 63 keystream bits which, assuming that the key has size 128, breaks the cipher.

Conclusion

In this paper we have described practical applications of certain cyclic codes over F q and F q n generated by a Boolean function. Determining the immunity against an algebraic attack on a WG-type construction involves determining the minimum value T s (n) where s are codewords of the cyclic codes generated by g 0 (x) = gcd(f (x), x q + x) and g 1 (x) = gcd(f (x) + 1, x q + x) over F q . This is a slightly di↵erent problem than finding minimum distances of a code and an algorithm for determining this ordered minimum distance problem is left as an open problem. Moreover, it does not seem that maximal algebraic immunity alone is su cient to ensure optimal values for T s (n). We propose that a strong Boolean function should attain maximal value among the codewords s in its two cyclic codes and that the analysis of this paper must be accounted for when designing secure word-based stream ciphers.

It should also be noted that our analysis may have applications to similar word-based constructions, most notably SNOW-3G [START_REF]Specification of the 3GPP Confidentiality and Integrity Algorithms UEA& UIA2 Document 2: Snow 3G Specification (version 1.1)[END_REF].

  the exponents in (3) must correspond to distinct integers of hamming-weight d. In a general sum of powers of L,f (L(x)) = P 2 k 1 i=0 L(x) i , the terms of weight d (3) can be reached by varying over each possible k d

  root of the LFSR feedback polynomial. Since the shift-register obeys a linear recursion, each bit s L(x↵ t+i ) where x is the initial state. If the state of the LFSR at time t is S

t+i of the state S t at time t can be described linearly by L t+i (x) = t = (s t , s t+1 , . . . , s t+n 1