Keywords: fonction non-lisse, méthode de Newton, non-convergence, M-matrice, P-matrice, problème de complémentarité linéaire

Non convergence de l'algorithme de Newton-min simple pour les problèmes de complémentarité linéaires avec P-matrice -Le rapport complet

Résumé : L'algorithme Newton-min, utilisé pour résoudre le problème de complémentarité linéaire (PCL) 0

x ⊥ (M x + q) 0 peut être interprété comme un algorithme de Newton non lisse sans globalisation cherchant à résoudre le système d'équations linéaires par morceaux min(x, M x + q) = 0, qui est équivalent au PCL. Lorsque M est une M-matrice d'ordre n, on sait que l'algorithme converge en au plus n itérations. Nous montrons dans cet article que ce résultat ne tient plus lorsque M est une P-matrice d'ordre n 3; l'algorithme peut en effet cycler dans ce cas. On a toutefois la convergence de l'algorithme pour une P-matrice d'ordre 1 ou 2.

Introduction

The linear complementarity problem (LCP) consists in finding a vector x 0 with n components such that M x + q 0 and x ⊤ (M x + q) = 0. Here M is a real matrix of order n, q is a vector in R n , the inequalities have to be understood componentwise, and the sign ⊤ denotes matrix transposition. The LCP is often written in compact form as follows LC(M, q) : 0 x ⊥ (M x + q) 0.

This problem is known to have a unique solution for any q ∈ R n if and only if M is a P-matrix [34,11], i.e., a matrix with positive principal minors: det M II > 0 for all nonempty I ⊂ {1, . . . , n}.

We denote by P the set of P-matrices. Other classes of matrices M intervening in the discussion below are the class Z of Z-matrices (which have nonpositive off-diagonal elements: M ij 0 for all i = j) and the class M := P ∩ Z of M-matrices (they are called K-matrices in [11]). Since the components of x and M x+q must be nonnegative in LC(M, q), their perpendicularity with respect to the Euclidean scalar product required in the problem is equivalent to the nullity of the Hadamard product of the two vectors, that is x • (M x + q) = 0.

(1.1)

Recall that the Hadamard product u • v of two vectors u and v is the vector having its ith component equal to u i v i . A point x such that (1.1) holds is here called a node or is said to satisfy complementarity. Since, for a node x, either x i or (M x + q) i vanishes, for all indices i, there are at most 2 n nodes for a nondegenerate matrix M , which is a matrix having nonsingular principal submatrices. On the other hand, a point x such that x and M x + q are nonnegative (resp. positive) is said to be feasible (resp. strictly feasible). A solution to LC(M, q) is therefore a feasible node. Many algorithms have been proposed to solve problem LC(M, q) [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]11]. They may be based on pivoting techniques [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF][START_REF] Cottle | Complementarity pivot theory of mathematical programming[END_REF], which often suffer from the combinatorial aspect of the problem (i.e., the 2 n possibilities to realize (1.1)), on interior point methods, which originate from an algorithm introduced by Karmarkar in linear optimization [21; 1984] (see also [23; 1991] for one of the first accounts on the use of interior point methods to solve linear complementarity problems), and on nonsmooth Newton approaches [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems (two volumes)[END_REF], such as the one considered here. See [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]11] for other iterative methods.

The algorithm we consider in this paper maintains the complementarity condition (1.1) at each iteration, while feasibility is obtained at convergence (when this one occurs). As a result, all the iterates are nodes, except possibly the first one, and the algorithm terminates as soon as it has found a feasible iterate. More specifically, suppose that the current iterate x is a node. The algorithm first defines index sets A + and I + associated with the next iterate x + ; in its simplest form, it takes A + := {i : x i (M x + q) i } and I + := {i : x i > (M x + q) i }.

(1.2)

Then it computes x + by solving the linear system formed of the equations

x + A + = 0 and (M x + + q) I + = 0. (1.3)
To have a well defined algorithm, an assumption on M is necessary so that this system has a solution for any choice of complementary sets A + and I + , and vector q: M must be nondegenerate. The motivation sustaining the algorithm is that it can be viewed as a semismooth Newton method to solve the system of piecewise linear equations min(x, M x + q) = 0, (1.4) in which the minimum operator 'min' acts componentwise: [min(x, y)] i = min(x i , y i ). On the one hand, since, for a and b ∈ R, min(a, b) = 0 if and only if 0 a ⊥ b 0, the system (1.4) is indeed equivalent to problem LC(M, q) (see [25; 1976] and [28; 1977, lemma 2.1] for instance). On the other hand, it is indeed clear that the function in (1.4) is differentiable at a point x without doubly active index (i.e., without index i such that x i = (M x + q) i ) and that its Jacobian matrix is the one used in the linear system (1.3); when there are doubly active indices, the Jacobian used in (1.3), determined by the choice (1.2), is an element of the Clarke generalized Jacobian [8] of the function in (1.4), so that the algorithm may indeed be viewed as a semismooth Newton method [33; 1993]. This description makes it natural to call Newton-min the algorithm that updates x by the formulas (1.2)- (1.3).

Here are some remarks about algorithm (1.2)-(1.3). First, note that the algorithm has a principle quite different from the one used by an interior point approach, which generates strictly feasible points, while complementarity is obtained at the limit. Note also that if, in the local analysis of the method, it is important to allow the first iterate x 1 not to be a node, in this paper, x 1 will always be assumed to be a node. Finally, observe that a consequence of the fact that the algorithm only generates nodes is that it is equivalent to say that it converges or that it converges in a finite number of iterations or that it does not cycle (algorithm (1.2)-(1.3) is a Markov process).

The algorithm sketched above and that we further explore in this paper can be traced back at least to the algorithm (6.2)-(6.4) of Aganagić in [1; 1984], who considers a linear complementarity problem that reads instead 0 (Xx) ⊥ (Y x+q) 0, in which X and Y must be jointly M-regular. Here, it is not important to be precise about the meaning of this property of joint M-regularity, but to know that if X = I, as in LC(M, q), then the property is equivalent to the M-matricity of Y ≡ M . When X = I and Y ≡ M , the algorithm (6.2)-(6.4) of Aganagić [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF] is exactly the algorithm (1.2)-(1.3) above, which is proven in [1; theorem 6.2] to be monotonically (in the sense that x k

x k+1 for k 2) and globally (i.e., x 1 may be arbitrary) convergent to the unique solution of LC(M, q), provided M ∈ M (in [START_REF] Aganagić | Newton's method for linear complementarity problems[END_REF], the first iterate x 1 is supposed to be zero, but this assumption is not necessary). Under the conditions that x 1 = 0 and M ∈ M, this algorithm is actually identical to the one proposed and analyzed earlier by Chandrasekaran [7; 1970] [1; remark 2]. Although not presented in that manner, the algorithm proposed by Bergounioux, Ito, and Kunisch [4;[START_REF] Bergounioux | Primal-dual strategy for constrained optimal control problems[END_REF] to solve a strongly convex quadratic optimal control problem, under the name of primal-dual active set strategy (PDAS), can be viewed as an extension of algorithm (1.2)-(1.3) to an infinite dimensional problem (see [3; 2000] for a comparison with interior point methods); the authors prove the convergence of the algorithm on a discretized version of the problem under some conditions. In [18; 2003], Hintermüller, Ito, and Kunisch consider a quadratic optimization problem with simple bounds, which can actually be written as a linear complementarity problem like LC(M, q); they establish the equivalence between the PDAS strategy and the semismooth Newton method, i.e., algorithm (1.2)-(1.3); the algorithm is also shown to be locally superlinearly convergent when M ∈ P, and to be monotonically (in the sense that

i x k i i x k+1 i
) and globally convergent when M is in the set that we denote here by M ε := {M : M is a matrix near an M-matrix of the same order}.

(1.5)

In M ε , the level of proximity to an M-matrix is left unprecise (and depends on the considered matrix M ), but in the proof of [18; theorem 3.4], this proximity is sufficiently tight to imply M ε ⊂ P. Another interesting property of algorithm (1.2)-(1.3), proved by Kanzow [20; 2004], is its convergence in at most n iterations when M ∈ M and the first iterate is a node. When applied to algorithm (1.2)-(1.3), the result of Fischer and Kanzow [15;1996] shows that a solution to LC(M, q), if any, is reached in one step, provided x 1 is sufficiently close to that solution and M is a nondegenerate matrix. To conclude this review of results, we would like to cite the quadratic local convergence of Newton's method for piecewise C 1 functions proved by Kojima and Shindo [24; 1986], which is related to the formulation (1.4) of problem LC(M, q). This paper presents examples of nonconvergence of the Newton-min algorithm when M is a P-matrix. These counter-examples hold for the plain (or undamped) Newton-min algorithm, i.e., without the use of globalization techniques such as linesearch or trust regions. One may believe INRIA that, as a Newton-like method to solve a nonlinear (and nonsmooth) system of equations, this is not a good strategy. We share this opinion, in general. However the algorithm deals with a piecewise linear function, and it has been shown to be convergent without globalization techniques when M is a P-matrix sufficiently near an M-matrix (see the discussion in the previous paragraph). Therefore, searching for the weakest assumptions for which these convergence properties hold seems to us a valid topic. The examples in this paper show that it is not enough to require the P-matricity for M .

The paper is structured as follows. In the next section, we are more specific on the definition of the algorithm, by making it a slightly more flexible than in the description (1.2)-(1.3) above. Some elementary properties of the algorithm, useful in the sequel, are also given. Section 3 describes and analyzes the examples of nonconvergence of the plain Newton-min algorithm with a P-matrix, when n 3. These counter-examples work for both definitions of the algorithm, those of sections 1 and 2. In them, the algorithm can be forced to cycle and visit p nodes, with a p that can be chosen arbitrarily in {3, . . . , n}. We consider successively the cases when n is odd and even, which require a different analysis. These counter-examples readily imply that the convergence radius of the plain Newton-min algorithm for solving LC(M, q) with M ∈ P, which is known to be > 0 [18; 2003], can be arbitrarily small (corollary 3.8). In section 4, the plain Newton-min algorithm is claimed to converge for a P-matrix when n = 1 or n = 2, a crumb of consolation. The paper concludes with the perspective section 5.

The algorithm

The plain Newton-min algorithm described in this section generates points that satisfy the complementarity conditions in (1.1), while the nonnegativity conditions in LC(M, q) are satisfied when the solution is reached. The starting point x 1 may or may not satisfy this complementarity condition.

Algorithm 2.1 (plain Newton-min) Let x 1 ∈ R n .
For k = 2, 3, . . ., do the following.

1. If x k-1 is a solution to LC(M, q), stop.

Choose complementary index sets A

k := A k 0 ∪ A k 1 and I k := I k 0 ∪ I k 1 ,
where

A k 0 := {i : x k-1 i < (M x k-1 + q) i }, A k 1 ⊂ E k := {i : x k-1 i = (M x k-1 + q) i }, I k 0 := {i : x k-1 i > (M x k-1 + q) i }, I k 1 := E k \ A k 1 .
3. Determine x k as a solution to

x k A k = 0 (M x k + q) I k = 0. (2.1)
The algorithm is well defined if M is nondegenerate. This assumption will be generally reinforced by supposing that M ∈ P. We recall from [START_REF] Fiedler | On matrices with nonpositive off-diagonal elements and principal minors[END_REF]11] that

M ∈ P ⇐⇒ any x verifying x • (M x) 0 vanishes. (2.2)
Note that algorithm 2.1 is more flexible than the one presented in section 1, in that the doubly active indices (those in E k ) can be chosen to belong either to A k or I k . If this choice is not entirely determined by the current point x k , the generated sequence {x k } may no longer be a Markov process. We will not be more precise here, however, since this choice has actually no impact on the counter-examples given in this paper, for which the algorithm never generates iterates with doubly active indices.

In this paper, we always assume that x 1 is a node; then there are two complementary subsets A 1 and I 1 of {1, . . . , n}, such that x 1 A 1 = 0 and (M x 1 + q) I 1 = 0. In other contexts, in particular for studying the local convergence of the algorithm [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF], it is better not to make this assumption.

By the selection of the index sets in step 2, one certainly has for k 2:

x k-1 A k (M x k-1 + q) A k and x k-1 I k (M x k-1 + q) I k . (2.3) Therefore, x k-1 A k 0 and (M x k-1 + q) I k 0 (2.4)
must hold [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Ch | Inexact semismooth Newton methods for large-scale complementarity problems[END_REF]. Indeed, for i

∈ A k , x k-1 i (M x k-1 + q) i by (2.
3) and, since either x k-1 i = 0 or (M x k-1 + q) i = 0 by complementarity, one necessarily has x k-1 i 0, which proves the first inequality in (2.4). A similar reasoning yields the second inequality in (2.4).

We denote by e i the ith vector of the canonical basis of R n : its jth component is equal to 1 if j = i and to 0 otherwise.

Nonconvergence for n 3

In this section, we show that the plain Newton-min algorithm described in section 2 may not converge if M is a P-matrix and n 3. We start in section 3.1 with the case when n is odd and provide an example of a P-matrix M , a vector q, and a starting point, for which the algorithm makes cycles having n nodes. In section 3.2, we consider the case when n is even and construct another class of P-matrices, which can be viewed as perturbations of the matrices in the first example, for which a cycle having n nodes is also possible. We conclude in section 3.3 by constructing examples for which the plain Newton-min algorithm makes cycles visiting p nodes, with p arbitrary in {3, . . . , n}.

Cycles with an odd number of nodes

In this section, we show the nonconvergence of the plain Newton-min algorithm for problems having the following features.

Example 3.1 Let n 2. The matrix M ∈ R n×n and the vector q ∈ R n are given by

M =       1 α α . . . . . . 1 α 1       and q = 1,
where 1 denotes the vector of all ones and the elements of M that are not represented are zeros.

More precisely,

M ij = 1 if i = j, M ij = α if i = (j mod n) + 1
, and M ij = 0 otherwise. Since q 0, a solution to LC(M, q) for these data is x = 0.

The matrix M and the vector q in example 3.1 have already been used by Morris [START_REF] Morris | Randomized pivot algorithms for P -matrix linear complementarity problems[END_REF] (in that paper, α = 2, M is the transpose of the one here, and q = -1), although we arrived at them in a INRIA different manner, as explained in section 4. For completeness and precision, we start by studying the P-matricity of the matrix M in example 3.1 (the result for n odd and α = 2 was already claimed in [START_REF] Morris | Randomized pivot algorithms for P -matrix linear complementarity problems[END_REF] without a detailed proof).

Lemma 3.2 Consider the matrix M in example 3.1. If n is even, M ∈ P if and only if |α| < 1. If n is odd, M ∈ P if and only if α > -1. Proof. Observe first that if α -1, the nonzero vector x = 1 is such that x•(M x) = (1+α)1 0; hence M / ∈ P. Observe next that M ∈ P when -1 < α < 0. Indeed, let x ∈ R n be such that x • (M x) 0 or equivalently x 1 (x 1 + αx n ) 0, x 2 (x 2 + αx 1 ) 0, . . . , x n (x n + αx n-1 ) 0. (3.1)
If x n > 0, using (3.1) from right to left shows that all the components of x are positive and verify

0 < x n |α|x n-1 |α| 2 x n-2 • • • |α| n-1 x 1 |α| n x n .
Since |α| < 1, this is incompatible with x n > 0. Having x n < 0 is not possible either (just multiply x by -1 and use the same argument). Hence x n = 0 and using (3.1) from left to right shows that x = 0. The P-matricity now follows from (2.2). When α = 0, M is the identity matrix and is therefore a P-matrix. Suppose now that n is even. If α 1, the nonzero vector x, defined by

x i = (-1) i+1 , is such that x • (M x) = (1 -α)1 0; hence M / ∈ P.
Let us now show, by contradiction, that M ∈ P if 0 < α < 1, assuming that there is a nonzero x satisfying x • (M x) 0, hence that (3.1) holds. Then, as above, all the components of x are nonzero and one can assume that x n > 0. Starting with the rightmost inequality of (3.1), one obtains by induction for i = 1, . . . , n -1:

x n-i -1 α i x n < 0 (for i odd) and

x n-i 1 α i x n > 0 (for i even).

Since n is even, x 1 -(1/α n-1 )x n < 0 and, using the first inequality of (3.1),

x n -(1/α)x 1 (1/α n )
x n , which is in contradiction with 0 < α < 1 and x n > 0. Suppose finally that n is odd and α > 0. Again, for proving that M ∈ P, we argue by contradiction, assuming that there is a nonzero x such that x • (M x) 0, which is equivalent to (3.1). As above, one can assume that x n > 0. Starting with the rightmost inequality in (3.1), one can specify by induction the sign σ(x i ) of the x i 's:

σ(x n-1 ) = -1, σ(x n-2 ) = 1, . . . , σ(x 1 ) = (-1) n-1 = 1,
since n is odd. Finally, the first inequality in (3.1) gives σ(x n ) = -σ(x 1 ) = -1, in contradiction with x n > 0. Hence x = 0 and M ∈ P by (2.2).

Recall that e k denotes the kth basis vector of R n . Lemma 3.3 Suppose that n 2 and consider problem LC(M, q) in which M and q are given in example 3.1 with α > 1. When applied to that problem LC(M, q) and started from x 1 = -e 1 , algorithm 2.1 cycles by visiting in order the n nodes x k = -e k , k = 1, . . . , n.

Proof. The proof proceeds by induction. By assumption, x 1 = -e 1 .

Suppose now that x k-1 = -e k-1 for some k ∈ {2, . . . , n} and let us show that x k = -e k . There hold

x k-1 =               0 . . . 0 -1 0 0 . . . 0               and M x k-1 + q =               1 . . . 1 0 1 -α 1 . . . 1              
,

where the component -1 (resp. 0) is at position k -1 in x k-1 (resp. in M x k-1 + q).
The update rules of algorithm 2.1 and α > 1 then imply that I k = {k} and x k = -e k . We still have to show that x n+1 = -e 1 . From x n = -e n , one deduces that

x n =        0 0 . . . 0 -1        and M x n + q =        1 -α 1 . . . 1 0        .
Again, the update rules of algorithm 2.1 and α > 1 show that I n+1 = {1} and x n+1 = -e 1 . Therefore algorithm 2.1 cycles.

By combining lemmas 3.2 and 3.3, one can see that algorithm 2.1 may cycle when M is a Pmatrix of odd order n 3: this is the case when the algorithm is started at x 1 = -e 1 and when M and q are given by example 3.1, with n odd and α > 1. When n is even, the condition α > 1 used in lemma 3.3 prevents the matrix M from being a P-matrix. It is not difficult to show, however, that the algorithm can also cycle when n is even, n 4, and M ∈ P, by using the construction given in the proof of proposition 3.7: the cycle visits an odd number of nodes (hence < n).

Cycles with an even number of nodes

The next family of examples is obtained by perturbing the matrices of example 3.1 with a parameter β, in order to construct cycles having n nodes for an even oder P-matrix. Example 3.4 Let n 3. The matrix M ∈ R n×n and the vector q ∈ R n are given by

M =          1 β α α . . . β β . . . 1 . . . α 1 β α 1          and q = 1,
where the elements of M that are not represented are zeros or, more precisely, for i, j ∈ {1, . . . , n}:

M ij =        1 if i = j α if i = (j mod n) + 1 β if i = ((j + 1) mod n) + 1 0 otherwise.
Since q 0, a solution to LC(M, q) for these data is x = 0.

INRIA

Conditions for having an n-node cycle with algorithm 2.1 on problem LC(M, q) with M and q given by example 3.4 are very simple to express. Lemma 3.5 Suppose that n 3 and consider problem LC(M, q) in which M and q are given in example 3.4 with α > 1 and β < 1. When applied to that problem LC(M, q) and started from x 1 = -e 1 , algorithm 2.1 cycles by visiting in order the n nodes defined by x k = -e k , k = 1, . . . , n.

Proof. The proof is quite similar to the one of lemma 3.3 and proceeds by induction. By assumption, x 1 = -e 1 .

Suppose now that x k-1 = -e k-1 for some k ∈ {2, . . . , n -1} and let us show that x k = -e k . There hold

x k-1 =                 0 . . . 0 -1 0 0 0 . . . 0                 and M x k-1 + q =                 1 . . . 1 0 1 -α 1 -β 1 . . . 1                
, where the component -1 (resp. 0) is at position k -1 in x k-1 (resp. in M x k-1 + q). The update rules of algorithm 2.1, α > 1, and β < 1 then imply that I k = {k} and x k = -e k . Therefore, by induction

x n-1 =          0 0 . . . 0 -1 0          and M x n-1 + q =          1 -β 1 . . . 1 0 1 -α          .
The update rules of algorithm 2.1, α > 1, and β < 1 then imply that I n = {n}, so that

x n =        0 0 . . . 0 -1        and M x n + q =          1 -α 1 -β 1 . . . 1 0          .
Again, the update rules of algorithm 2.1, α > 1, and β < 1 show that I n+1 = {1}. Therefore, x n+1 = -e 1 and algorithm 2.1 cycles.

We have now to examine whether the conditions on α and β given by lemma 3.5 are compatible with the P-matricity of the matrix M in example 3.4. Actually, the conditions on α and β ensuring the P-matricity of that matrix M are nonlinear and much more complex to write than the simple conditions on α given in lemma 3.2; in particular, their number depends on the dimension n. The shaded regions in figure 3.1 show the intersections with the box [-2, 2] × [-2, 2] of the sets P formed of the (α, β) pairs for which the matrix M is a P-matrix, when its order is n = 3, 4, 5, and 6. Of course, by lemma 3.2, P contains the set {(α, β) : -1 < α < 1, β = 0} when n is even and the set {(α, β) : -1 < α, β = 0} when n is odd. It is not clear at this point, however, whether, for any n 3, these regions will contain points with α > 1 and β < 1, which are the conditions highlighted by lemma 3.5. Lemma 3.6 below shows that this is actually the case since P always contains the interiors of the nonconvex polyhedron and the ellipse represented in figure 3.1, which are independent of n. 2] for which the matrix M of example 3.4 is a P-matrix, when its order is n = 3, 4, 5, and 6; these regions are nonconvex but star-shaped with respect to (0, 0), which is the point corresponding to the identity matrix. According to lemma 3.6, the interiors of the represented nonconvex polyhedron and ellipse, which are independent of n, are always contained in these regions, for any n 3.

To prepare the proof of lemma 3.6, we write the circulant matrix M in example 3.4 as follows

M = I + βJ n-2 + αJ n-1 ,
where I denotes the n × n identity matrix and J is the elementary circulant n × n matrix

J =       0 1 0 0 . . . . . . 1 1 0      
.

INRIA

More precisely, J ij = 1 if j = (i mod n) + 1 and J ij = 0 otherwise. It is well known [16; formula (4.7.10)] that J is diagonalizable on C, meaning that there is a diagonal matrix D ∈ C n×n and a nonsingular matrix P ∈ C n×n such that J = P DP -1 ; in addition its eigenvalues are the nth roots of unity: D := Diag(1, w 1 , . . . , w n-1 ), where w = e 2πi/n and i is the imaginary unit.

Lemma 3.6 The set of (α, β) pairs ensuring that M ∈ P in example 3.4 contains the set

{(α, β) : |α| -1 < β < |α|/4 or α 2 + 8(β -1 2 ) 2 < 2}.
Proof. We only have to show that when α and β satisfy

|α| -1 < β < |α| 4 or α 2 + 8 β - 1 2 2 < 2, ( 3.2)
the matrix M + M ⊤ is positive definite, since then M is clearly a P-matrix [22; p. 175].

For any integer p, (J p ) ⊤ = J n-p . Therefore M + M ⊤ = 2I + αJ + βJ 2 + βJ n-2 + αJ n-1 and, using J = P DP -1 , we obtain

M + M ⊤ = P 2I + αD + βD 2 + βD n-2 + αD n-1 P -1 .
This identity shows that the eigenvalues of the symmetric matrix M + M ⊤ are the (necessarily real) numbers λ k := 2 + αe 2kπi/n + βe 4kπi/n + βe 2k(n-2)πi/n + αe 2k(n-1)πi/n , for k = 0, . . . , n -1. Using e 2kpπi/n + e 2k(n-p)πi/n = 2 cos(2kpπ/n) (for p integer) and cos 2θ = 2 cos 2 θ -1, we obtain

λ k = 2 + 2α cos(2kπ/n) + 4β cos 2 (4kπ/n) -2β.
We see that the desired positivity of the eigenvalues λ k depends on the positivity of the following polynomial on [-1, 1]:

t → ϕ(t) = 2βt 2 + αt + (1-β).
We denote by t

± := [-α ± (α 2 -8β(1 -β)) 1/2
]/(4β) the roots of ϕ when β = 0 and consider in sequence the three possible cases, identifying in each case conditions that ensure the positivity of ϕ on [-1, 1].

• Case β = 0. Then ϕ is positive on [-1, 1] if |α| < 1.
• Case β > 0. There are two subcases.

• If ϕ has no real root, i.e., if

α 2 < 8β(1 -β) or equivalently α 2 + 8(β -1 2 ) 2 < 2, then ϕ is clearly positive on R. • If ϕ has two (possibly equal) real roots t ± , i.e., if α 2 8β(1 -β), then these verify t -t + and ϕ is positive on [-1, 1], -either if 1 < t -, which is ensured if -α -1 < β < -α/4, -or if t + < -1, which is ensured if α -1 < β < α/4.
• Case β < 0. Then, ϕ has two real roots t ± , which verify t + < t -, and ϕ is positive on [-1, 1] if both t + < -1 and 1 < t -holds, which is ensured if |α| -1 < β < 0.

By gathering the above conditions, we obtain (3.2).

RR n°7160

Nonconvergence

The nonconvergence result is summarized in proposition 3.7, in which it is shown that cycles made of p nodes are possible when n 3 and p ∈ {3, . . . , n}. It is then shown with proposition 3.8 that when n 3, although the plain Newton-min algorithm is known to converge locally, i.e., when the starting point is near the solution to LC(M, q), the radius of convergence can be made as small as desired by modifying q; this makes the convergence of the plain Newton-min algorithm unlikely.

Proposition 3.7 (nonconvergence for n 3) When n 3, algorithm 2.1 may fail to converge when trying to solve LC(M, q) with a P-matrix M . A cycle made of p nodes is possible, for an arbitrary p ∈ {3, . . . , n}.

Proof. Since the plain Newton-min algorithm visits only a finite number of nodes, it fails to converge if and only if it cycles.

When n 3 is odd, a cycle made of n nodes is possible on problem LC(M, q) with M and q given by example 3.1, and α > 1: lemma 3.2 shows that M ∈ P and lemma 3.3 shows that a cycle is possible.

When n 4 is even, a cycle made of n nodes is possible on problem LC(M, q) with M and q given by example 3.4, and α and β satisfying α > 1 and α -1 < β < α/4 (hence β < 1/3): lemma 3.5 shows that a cycle is possible and lemma 3.6 shows that M ∈ P.

When n 3 and p ∈ {3, . . . , n}, consider a p × p matrix M ∈ P, a vector q ∈ R p , and a starting point x1 ∈ R p , such that algorithm 2.1 applied to problem LC( M , q) and starting at x1 generates iterates xk forming a cycle made of p nodes (this is possible by what has just been proven). With obvious notation, define

M = M 0 p×(n-p) 0 (n-p)×p I n-p , q = q 0 n-p
, and x 1 = x1 0 n-p .

The P-matricity of M is clear, by observing that x • (M x) 0 implies x = 0. Denote by x k the kth iterate generated by algorithm 2.1 on LC(M, q) starting from x 1 . Observe first that when an index i > p, there holds x k i = 0, whenever i ∈ I k or i ∈ A k ; therefore the generated iterates x k ∈ R p × {0 n-p }. Hence, if the same rule as the one used by algorithm 2.1 on R p is used to decide whether an index i ∈ E k will be considered as being in I k or A k , the iterates x k will be (x k , 0 n-p ). Obviously, as the xk 's, these iterates also form a cycle made of p nodes.

To make the above nonconvergence result more concrete, we provide two examples of Pmatrices M n , of order n = 3 and n = 4 respectively, which make algorithm 2.1 fail with an n-cycle, when it starts at x 1 = (-1, 0, . . . , 0) to solve problem LC(M n , 1):

M 3 :=   1 0 2 2 1 0 0 2 1   and M 4 :=     1 0 1/2 4/3 4/3 1 0 1/2 1/2 4/3 1 0 0 1/2 4/3 1     . (3.3)
We have used the lemmas 3.2 and 3.3 for constructing M 3 and the lemmas 3.5 and 3.6 for designing M 4 .

The convergence radius of an iterative algorithm for finding a solution x to a given problem is the largest ρ > 0 such that the algorithm converges to x when the initial iterate is in the ball of center x and radius ρ. The plain Newton-min algorithm is known to be locally convergent [18; 2003]. By scaling the variables in the previous examples, however, one can make the convergence radius of the plain Newton-min algorithm as small as desired. Even though this observation is straightforward, INRIA Proposition 4.2 (convergence for n = 1) Suppose that n = 1, that M is nondegenerate (M = 0), and that LC(M, q) has a solution. Then the plain Newton-min algorithm converges.

Proof. Without restriction, it can be assumed that the first iterate x 1 is a node. If x 1 is a solution, the algorithm stops at that point (by step 1 of algorithm 2.1). If x 1 is not the solution, there is another node that is solution (by assumption). Since there are no more than 2 nodes (since n = 1), the algorithm takes the solution as the next iterate (by lemma 4.1) and stops there.

Since a P-matrix is nondegenerate, the elementary preceding result applies when M ∈ P. When M = 0 ∈ R but LC(M, q) has no solution, the algorithm cycles by visiting alternatively the 2 nodes of the problem, i.e., x = -q/M < 0 and x = 0.

We start the study of the convergence of the plain Newton-min algorithm when n = 2 by showing that the algorithm cannot do cycles made of 2 nodes (lemma 4.3). We have seen with examples 3.1 and 3.4 that the algorithm can do a 3-cycle when n

3, but this implies some conditions that are highlighted by lemma 4.4. We finally show that these conditions cannot be satisfied when n = 2 and, as a result, that the algorithm must converge (proposition 4.5).

Lemma 4.3 (no 2-cycle) If M is a P-matrix, then the plain Newton-min algorithm does not make cycles formed of 2 distinct nodes.

Proof. We argue by contradiction, assuming that the algorithm visits in order the following nodes x 1 → x 2 → x 1 , with x 1 = x 2 . Since the algorithm goes from x 1 to x 2 and from x 2 to x 1 , the definition of the sets A k and I k in step 2 of the algorithm implies that

x 1 A 2 (M x 1 + q) A 2 and x 1 I 2 (M x 1 + q) I 2 , (4.1) x 2 A 1 (M x 2 + q) A 1 and x 2 I 1 (M x 2 + q) I 1 . (4.2)
After possible rearrangement of the component order, we get

x 2 -x 1 =     0 A 1 ∩A 2 x 2 A 1 ∩I 2 0 I 1 ∩A 2 x 2 I 1 ∩I 2     -     0 A 1 ∩A 2 0 A 1 ∩I 2 x 1 I 1 ∩A 2 x 1 I 1 ∩I 2     =     0 A 1 ∩A 2 x 2 A 1 ∩I 2 -x 1 I 1 ∩A 2 (x 2 -x 1 ) I 1 ∩I 2     .
Observe that the components of

x 2 -x 1 with indices in A 1 ∩ I 2 are nonpositive since x 2 A 1 ∩I 2 (M x 2 + q) A 1 ∩I 2 [by (4.2) 1 ] = 0 [by (2.1) 2 ]
and that the components with indices in I 1 ∩ A 2 are nonnegative since -x 1

I 1 ∩A 2 -(M x 1 + q) I 1 ∩A 2 [by (4.1) 1 ] = 0 [by (2.1) 2 ]
. On the other hand, by (2.1) 2 , there holds

M (x 2 -x 1 ) =     (M x 2 ) A 1 ∩A 2 -q A 1 ∩I 2 (M x 2 ) I 1 ∩A 2 -q I 1 ∩I 2     -     (M x 1 ) A 1 ∩A 2 (M x 1 ) A 1 ∩I 2 -q I 1 ∩A 2 -q I 1 ∩I 2     =     (M (x 2 -x 1 )) A 1 ∩A 2 -(M x 1 + q) A 1 ∩I 2 (M x 2 + q) I 1 ∩A 2 0     .
In this vector, the components with indices in

A 1 ∩ I 2 are nonnegative since -(M x 1 + q) A 1 ∩I 2 -x 1 A 1 ∩I 2 [by (4.1) 2 ] = 0 [by (2.1) 1 ] and the components with indices in I 1 ∩ A 2 are nonpositive since (M x 2 + q) I 1 ∩A 2 x 2 I 1 ∩A 2 [by (4.2) 2 ] = 0 [by (2.1) 1 ]. Therefore (x 2 -x 1 ) • M (x 2 -x 1 ) 0.
INRIA Since M ∈ P, there holds x 1 = x 2 by (2.2), contradicting the initial assumption.

Lemma 4.4 (necessary conditions for a 3-cycle) Suppose that M is a P-matrix and that the plain Newton-min algorithm cycles by visiting in order the three distinct nodes x 1 → x 2 → x 3 . Then the following three sets of indices must be nonempty

   (A 1 ∩ I 2 ∩ I 3 ) ∪ (I 1 ∩ A 2 ∩ A 3 ) (I 1 ∩ A 2 ∩ I 3 ) ∪ (A 1 ∩ I 2 ∩ A 3 ) (I 1 ∩ I 2 ∩ A 3 ) ∪ (A 1 ∩ A 2 ∩ I 3 ). (4.3)
Proof. We use the same technique as in the proof of lemma 4.3. Since the algorithm goes from x 1 to x 2 and from x 2 to x 3 , one has from step 2 of the algorithm that

x 1 A 2 (M x 1 + q) A 2 and x 1 I 2 (M x 1 + q) I 2 (4.4) x 2 A 3 (M x 2 + q) A 3 and x 2 I 3 (M x 2 + q) I 3 . (4.5)
Using (2.1) 1 , there holds

x 2 -x 1 =              0 A 1 ∩A 2 ∩A 3 0 A 1 ∩A 2 ∩I 3 x 2 A 1 ∩I 2 ∩A 3 x 2 A 1 ∩I 2 ∩I 3 0 I 1 ∩A 2 ∩A 3 0 I 1 ∩A 2 ∩I 3 x 2      . [+] [+] [+] [-] [0] [0]
The sign of the components given in the extra column on the right is justified as follows:

(M x 1 + q) A 1 ∩I 2 x 1 A 1 ∩I 2 = 0 [by (4.4) 2 and (2.1) 1 ] (M x 2 + q) I 1 ∩A 2 ∩A 3 x 2 I 1 ∩A 2 ∩A 3 = 0 [by (4.5) 1 and (2.1) 1 ] (M x 2 + q) I 1 ∩A 2 ∩I 3 x 2 I 1 ∩A 2 ∩I 3 = 0 [by (4.5) 2 and (2.1) 1 ].
Taking the Hadamard product of the two vectors now gives

(x 2 -x 1 ) • M (x 2 -x 1 ) =              0 A 1 ∩A 2 ∩A 3 0 A 1 ∩A 2 ∩I 3 -x 2 A 1 ∩I 2 ∩A 3 • (M x 1 + q) A 1 ∩I 2 ∩A 3 -x 2 A 1 ∩I 2 ∩I 3 • (M x 1 + q) A 1 ∩I 2 ∩I 3 -x 1 I 1 ∩A 2 ∩A 3 • (M x 2 + q) I 1 ∩A 2 ∩A 3 -x 1 I 1 ∩A 2 ∩I 3 • (M x 2 + q) I 1 ∩A 2 ∩I 3 0 I 1 ∩I 2 ∩A 3 0 I 1 ∩I 2 ∩I 3              , [0] [0] [-] [+] [+] [-] [0] [0]
where the extra column on the right gives the sign of each components. Therefore, if the index set (A 1 ∩ I 2 ∩ I 3 ) ∪ (I 1 ∩ A 2 ∩ A 3 ) were empty, we would have (x 2 -x 1 ) • M (x 2 -x 1 ) 0, which, with the P-matricity of M and (2.2), would imply that x 2 = x 1 , in contradiction with the initial assumption. We have proven that the first index set in (4.3) is nonempty.

To show that the second index set in (4.3) is nonempty, we use the fact that the algorithm goes from x 2 to x 3 and from x 3 to x 1 . Therefore, by cycling the indices in the result just obtained, we see that (A 2 ∩ I 3 ∩ I 1 ) ∪ (I 2 ∩ A 3 ∩ A 1 ) = ∅; this corresponds to the second index set in (4.3). By cycling the indices again, we obtain (A 3 ∩ I 1 ∩ I 2 ) ∪ (I 3 ∩ A 1 ∩ A 2 ) = ∅; this corresponds to the third index set in (4.3). Proposition 4.5 (convergence for n = 2) Suppose that M is a P-matrix and that n = 2.

Then the plain Newton-min algorithm converges.

Proof. We know that the algorithm converges if it does not make a p-cycle, a cycle made of p 2 distinct nodes. It cannot make a 2-cycle by lemma 4.3. By lemma 4.4, to make a 3-cycle, the three sets in (4.3) must be nonempty; but these sets are disjoint; since n = 2, one of them must be empty, so that the algorithm does not make a 3-cycle. Therefore, the algorithm will visit a 4th node if it has not found the solution on the first 3 visited nodes. This last node is then the solution, since the solution exist (M ∈ P) and there are at most 2 n = 4 distinct nodes (M ∈ P).

Perspectives

The work presented in this paper can be pursued along at least two directions. One possibility is to better mark out the set of matrices for which the plain Newton-min method converges. According to [18; 2003, theorem 3.4] and the counter-examples of the present paper, this set contains the set M ε of matrices sufficiently near an M-matrix but not all of the larger set P of P-matrices, see the discussion around (1.5). Such a study may result in the identification of an analytically well defined set of matrices or it may be a long process with an endless refinement. In the first case, it would be nice to see whether membership in that new matrix class can be determined in polynomial time, knowing that recognizing a P-matrix is a co-NP-complete problem [START_REF] Coxson | The P -matrix problem is co-NP-complete[END_REF][START_REF] Tseng | Co-NP-completeness of some matrix classification problems[END_REF]. Another possibility is to modify the algorithm to force its convergence for P-matrices or an even larger class of matrices. Being able to deal with P-matrices is important for at least three reasons.

INRIA

6 Figure 3 . 1 :

 631 Figure 3.1: The shaded regions are formed of the(α, β) pairs in [-2, 2] × [-2,2] for which the matrix M of example 3.4 is a P-matrix, when its order is n = 3, 4, 5, and 6; these regions are nonconvex but star-shaped with respect to (0, 0), which is the point corresponding to the identity matrix. According to lemma 3.6, the interiors of the represented nonconvex polyhedron and ellipse, which are independent of n, are always contained in these regions, for any n 3.

Example 3 . 1 ∩ A 2 ∩ A 3 =

 3123 1 was actually obtained for n = 3, by forcing the 3 sets in (4.3) to be nonempty by settingI {1}, A 1 ∩ I 2 ∩ A 3 = {2}, and A 1 ∩ A 2 ∩ I 3 = {3},which yields I k = {k}.

RR n°7160

Acknowledgments

We would like to thank P. Knabner and S. Kräutle for having drawn our attention on the Newtonmin algorithm and its efficiency on practical problems [6, 26; 2009, 2010], and M. Hintermüller, Ch. Kanzow, N. Metla [30; 2008], and A. Griewank for various exchanges of mails and conversations on the topics. Also, we wish to thank the two referees for their helpful comments. This work was partially supported by the MoMaS group (PACEN/CNRS, ANDRA, BRGM, CEA, EDF, IRSN).

we highlight it in the following corollary to stress the fact that without modification the plain Newton-min algorithm may have little chance to converge. Corollary 3.8 (small convergence radius for n 3) When n 3, the convergence radius of algorithm 2.1 to solve LC(M, q) with a P-matrix M may be arbitrarily small. Proof. Take an example of matrix M ∈ P and vector q 0 such that algorithm 2.1 cycles when it tries to solve LC(M, q) from some nonzero x 1 (this is possible by using one of the problems considered in the proof of proposition 3.7). The convergence radius of algorithm 2.1 for that problem is therefore less than x 1 (the solution is 0). Now, algorithm 2.1 starting at x1 = εx 1 , for some ε > 0, for solving LC(M, εq) generates the iterates xk = εx k and therefore cycles. The convergence radius for the plain Newton-min algorithm on this new problem is less than ε x 1 , which can be made arbitrarily small by letting ε ↓ 0.

4 Convergence for n = 1 or 2

In this section, we prove the convergence of the plain Newton-min algorithm, algorithm 2.1, when M is a P-matrix of order 1 or 2. The proof for n = 1 is straightforward. The one presented for n = 2 is indirect but highlights the origin of the counter-example 3.1 and allows us to present some properties of the plain Newton-min algorithm.

The convergence of the plain Newton-min algorithm when n = 1 is a direct consequence of the following reassuring and elementary property, already proven in [4; theorem 2.1] in the case where the complementarity problem expresses the optimality conditions of an infinite dimensional quadratic optimization problem.

Lemma 4.1 (stagnation at a solution) Suppose that M is nondegenerate. Then, a node is a solution to LC(M, q) if and only if the plain Newton-min algorithm, without its stopping test in step 1, starting at that node, takes the same node as the next iterate.

Proof. Let x be a node of problem LC(M, q). We denote by A + and I + the index sets determined in step 2 of algorithm 2.1 and by x + the next iterate.

If x is a solution, then x 0 and M x + q 0. There hold x + A + = 0 by (2.1) and x A + = 0 by (2.4) and the nonnegativity of x, so that x + A + = x A + . Similarly, there hold (M x + + q) I + = 0 by (2.1) and (M x + q) I + = 0 by (2.4) and the nonnegativity of M x + q, so that 0 = (M (x + -x)) I + = M I + I + (x + -x) I + [since (x + -x) A + = 0] and therefore (x + -x) I + = 0 by the nonsingularity of M I + I + . We have shown that x + = x.

Conversely, assume that x + = x. By (2.3), there hold

] and (M x + q) I + = (M x + + q) I + = 0 [by (2.1)], we get x 0 and M x + q 0; hence x is a solution.

We consider now the case when n = 1 and M is nondegenerate. For such a trivial problem, a direct proof of convergence is obviously possible. The proof given below uses instead lemma 4.1, which is nevertheless useful elsewhere.

First, linear complementarity problems with a P-matrix are encountered in practice [35; 2004]. Next, this is exactly the class of matrices that ensure the existence and uniqueness of the solution to the LCP [34,11], which forces us to pay attention to these matrices. Finally, the possibility to find a polynomial algorithm to solve the LCP with a P-matrix still seems to be an open question. Some authors argue that such an algorithm might exist; see Morris [31; 2002], who refers to a contribution by Megiddo [29;1988], himself citing an unpublished related note of Solow, Stone and Tovey [36; 1987]. The possibility that a modified version of the plain Newton-min algorithm might have the desired polynomiality property cannot be excluded. A natural remedy would be to add a globalization technique (linesearch or trust regions, see [START_REF] Conn | Trust-Region Methods[END_REF][START_REF] Bonnans | Numerical Optimization -Theoretical and Practical Aspects[END_REF] for example) to the plain Newton-min algorithm in order to force its convergence; see [17, 19; 1990, 2009], for contributions along that direction. This globalization is not straightforward since the Newton-min direction may not be a descent direction of the standard ℓ 2 merit function associated with (1.4) [2].


Claude Marché

Asma Tafat

Weakest Precondition Calculus, Revisited using Why3

Keywords: Deductive Verification, Weakest Precondition Calculus, Blocking Semantics, Automated Provers, Coq Proof Assistant déductive, calcul de plus faible pré-condition, sémantique bloquante, prouveurs automatiques, assistant | Ominus | Omult | Ole ( ** operators +, -, *, ≤

This report has two objectives. First, we present an original method of proof of soundness of a weakest precondition calculus, based on the notion of blocking semantics. The method mimics, at the level of logic specifications, the classical proof of type soundness. Moreover, the proof is performed formally using the Why3 environment for deductive verification, and we illustrate, along the development of the case study, the advanced features of Why3 we used. The result is a revisited presentation the weakest precondition calculus which is easy to follow, although formally made, thanks in particular to the high degree of proof automation that allows us to focus on the key points.

List of Figures

The Why3 system [START_REF] Bobot | Why3: Shepherd your herd of provers[END_REF] is an environment for deductive program verification. It allows the user to write programs annotated with logical specifications, and then generate proof obligations that guarantee the functional correctness.

Why3 can submit these obligations to a large collection of provers: automatic ones such as SMT solvers (Alt-Ergo [START_REF] Bobot | The Alt-Ergo automated theorem prover[END_REF], CVC3 [START_REF] Barrett | CVC3[END_REF], Z3 [START_REF] De Moura | Z3, an efficient SMT solver[END_REF], etc.) or interactive proof assistants like Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development[END_REF] or PVS [START_REF]The PVS system[END_REF].

In the spectrum of deductive verification tools, Why3 positions itself in the middle of interactive environments, that offer highly expressive languages but a low level of automation (Coq, PVS, Isabelle/HOL, etc.), and more automated systems but equipped with poorer languages (Spec#, VCC, Frama-C, KeY, etc.). A natural objective is to catch the best possible between these extremes: a sufficiently expressive language that still permits significantly automated proofs.

A first objective of this report is to illustrate the capabilities of Why3 on some case study, that involves complex objects like recursive algebraic data types and inductive predicate definitions, so as to illustrate how proofs can be automatized in such a context. A class of examples involving such kind of objects is given by the symbolic computation methods, in particular the tools for analyzing programs: compilers, static analyzers, etc. We chose to formalize a calculus of weakest precondition, partially following a recent approach due to Herms et al. [START_REF] Herms | A certified multi-prover verification condition generator[END_REF] for the certification of a verification condition generator, using Coq. This new approach is based on a notion of big-step blocking semantics, defined co-inductively. The second objective of this work is to present an original approach of the weakest precondition calculus. The originality resides in the use of a small-steps blocking semantics, which is moreover formalized using the Why3 language, significantly restricted with respect to those of Coq.

Apart from the algebraic data types and the inductive predicates, the advances features of Why3 that we are going to use concern the tools available to build proofs. First, we use a transformation of Why3 to make structural inductions. Second, when the goals to prove get too complex to be proved fully automatically, we call the Coq back-end, and within Coq we use the « why3 » tactic that can finish the proof of a sub-goal by calling automated provers through Why3. These features allow us to make in Coq only the subtle parts of the proofs.

Our case study is structured into components, named theories in Why3, and this report follows the same structure, represented on Figure 1. In Section 2 we define the language under study: its syntax (Section 2.1), its operational semantics (Section 2.2) and its typing rules (Section 2.3). In Section 2.4 we formulate and prove classical results of preservation of typing by reduction. Section 3 is devoted to the notion of substitution and to the important problems related to fresh variables. Although not needed for the remaining, Section 4 presents the classical Hoare logic rules and prove their soundness. The main results are then presented in Section 5: definition of the weakest precondition calculus, and soundness results based on the small-steps blocking semantics.

Formalization of Toy Imperative Language

This section formalize our language which contains assignment, sequence, a conditional statement and a while loop.

For each of the theories we will present, we are going to tell how the proof of the lemmas are obtained, by giving the table that is automatically generated by Why3 from the proof sessions [START_REF] Bobot | The Why3 platform, version 0.80[END_REF].

In most cases, provers are launched with a time limit of 5 seconds and a memory limit of 1 gigabyte. The cases where either one of these limits is reach are indicated in the tables below by results between parentheses.

Syntax

The first theory defines our language syntax. Traditionally, one defines language syntax with grammars. Here, we directly give Why3 definitions under the form of algebraic data type declarations. We can not present the Why3 syntax in detail in this report, you may refer to the manual [START_REF] Bobot | The Why3 platform, version 0.80[END_REF] or J.-C.Filliâtre lecture given at JFLA 2012 1 . The syntax is similar to the OCaml one, we hope it is intuitive enough for the reader.

Our language has three base types which correspond to three sets of values. We introduce abstract types for identifiers, by making a distinction between identifiers for mutable variables, and identifiers for logic variables.

Formulas of the logic language include terms (if Boolean), classical connectors for conjunction, negation and implication, and universal quantification. These are the only needed connectors, but others could be added. We also add local binding with let that allows us to write more elegant rules in the weakest precondition calculus. Note that we do not use any advanced technique such as De Bruijn indexes for local bindings. Difficulties related to variable naming will be handled in Section 3. We consider assignment of a mutable variable by a term, sequences (skip constant representing an empty sequence), standard conditional statement and a while loop.

Unconventionally, we add logic annotations directly in the syntax: explicit invariants for loops, as well as assert statements.

This first theory contains only definitions and no lemmas. Therefore, there is nothing to prove.

Operational semantics

The aim of the second theory is to define the operational semantics of our language. We choose a small step version.

In a first step, we define evaluation environments: for mutable variable, on one hand, and for logic variables on the other. get _ stack function provides a shortcut to access to an element; it is defined by induction on the list. Note that this function returns void if we try to access to a variable absent from the list. This case never happen on well-typed programs as defined in the next section. Finally, we pose lemmas showing a similar behavior to maps, that are automatically proved.

lemma get _ stack _ eq: forall x:ident, v:value, pi:stack. get _ stack x (Cons (x,v) pi) = v lemma get _ stack _ neq: forall x i:ident, v:value, pi:stack.

x = i → get _ stack i (Cons (x,v) pi) = get _ stack i pi
Next, we define an auxiliary function for evaluating binary operations. This definition is made by patternmatching on the values and the operator. Since in Why3 functions must be total, we complete our definition by returning void in ill-typed cases. 2 . We can now define evaluation of terms as a total recursive function that returns a value, and evaluation of formulas by a total recursive predicate. These two functions take a current state of a program, given by two environments Σ and Π, as arguments. → eval _ fmla sigma pi f1 ∧ eval _ fmla sigma pi f2

| Fnot f → not (eval _ fmla sigma pi f) | Fimplies f1 f2 → eval _ fmla sigma pi f1 → eval _ fmla sigma pi f2 | Flet x t f → eval _ fmla sigma (Cons (x,eval _ term sigma pi t) pi) f | Fforall x TYint f → forall n:int. eval _ fmla sigma (Cons (x,Vint n) pi) f | Fforall x TYbool f → forall b:bool. eval _ fmla sigma (Cons (x,Vbool b) pi) f | Fforall x TYunit f → eval _ fmla sigma (Cons (x,Vvoid) pi) f end
Note that in the case of binders Flet and Fforall, the linked variable x is added to the stack. In the following, we will note these functions in the form t Σ,Π and f Σ,Π .

Finally, we define the notion of valid formula, in the sense that it is true in any state.

predicate valid _ fmla (p:fmla) = forall sigma:env, pi:stack. eval _ fmla sigma pi p

Execution of a statement in a given state is classically defined by a one step reduction relation noted

Σ, Π, s Σ , Π , s
Since there are statements that do not reduce, we cannot formalize this relation by a Why3 function 3 . Instead, we use an inductive predicate of Why34 . Our definition is a formalization close to the usual representation using inference rules. It is shown on Figure 2. Thus, the one _ step _ while _ true case correspond to the rule inv Σ,Π cond Σ,Π = T rue Σ, Π, while cond invariant inv do body Σ, Π, body; while cond invariant inv do body

We use a blocking semantic in the way proposed by Herms et al. [START_REF] Herms | A certified multi-prover verification condition generator[END_REF]: The execution blocks whenever an invalid assertion is met or if a loop invariant is broken when the loop condition is checked (as shown in the rule above). Conversely, assuming that the program is well-typed, the only reason why a statement would not execute is that one of its annotations is not respected.

The advantage of such a definition is that we define the fact that a program respect its specifications by the fact that it executes without blocking [START_REF] Herms | A certified multi-prover verification condition generator[END_REF], and this applies in particular for a program does not terminate.

Another inductive predicate defines the execution in n steps of a program, by a reflexive-transitive closure, shown on Figure 3. Remark that we make explicit the number n of steps of execution in this definition, to do proofs by induction on n. Notice that the built-in int type in Why3 is the type of integers (signed), there is no built-in type for natural numbers only. Thus, to make well-founded inductions, we pose the following elementary, but essential, lemma. To end this theory, we define the notion of reducibility of a program in a given state.

predicate reducible (sigma:env) (pi:stack) (s:stmt) = exists sigma':env, pi':stack, s':stmt. one _ step sigma pi s sigma' pi' s'

So, for the case of our theory modeling the operational semantics, there are 3 lemmas to prove. Results are shown on Figure 4. As expected, the lemma steps _ non _ neg cannot be proved by automated provers, as it requires an induction on the many _ steps predicate. We do the proof in Coq, and this proof is easy (one line): induction 1; auto with zarith.

Typing

This theory introduces typing of programs. At first, we define a total function that returns a type of some value.

function type _ value (v:value) : datatype = match v with | Vvoid → TYunit | Vint _ → TYint | Vbool _ → TYbool end
Type of the binary operators of our language is given by a simple inductive predicate.

Relations between typing and operational semantic

This theory contains our first results linking typing and execution. At first, it defines a predicate for compatibility between typing and evaluation environments, that is, to each identifier of type t is associated a value of type t.

predicate compatible _ env (sigma:env) (sigmat:type _ env) (pi:stack) (pit: type _ stack) = (forall id: mident. type _ value (get _ env id sigma) = get _ reftype id sigmat) ∧ (forall id: ident. type _ value (get _ stack id pi) = get _ vartype id pit)

Then we pose an inversion lemma, that allows us to know which is the constructor of a value, from its type.

lemma type _ inversion : forall v:value. match (type _ value v) with | TYbool → exists b: bool. v = Vbool b | TYint → exists n: int. v = Vint n | TYunit → v = Vvoid end
The results of proofs of all the lemmas in this theory are given on Figure 8. For this particular lemma, an obvious case-based reasoning allows to prove such lemma. Nevertheless, theorem provers are not able to prove this goal 5 . We could do this proof with Coq, but we can also use Why3 features: goals transformations. In Why, there exists such a transformation to make structural induction. Here, the type value is not recursive, so it amounts to a case-based reasoning. The three generated sub-goals are automatically discharged. We also made a Coq proof of type _ inversion lemma. It can be done easily thanks to another feature of Why3: the why3 tactic of Coq. It is a tactic loaded as Coq plug-in by Require Import Why3. Ltac ae := why3 "alt-ergo" timelimit 5

Proof obligations

Alt-Ergo (0.93. The tactic ae is therefore a shortcut to call Alt-Ergo from Coq, with a timeout fixed to 5 seconds. The tactic allows, then, to call (through Why3) theorem provers, according to the architecture of Figure 9.

The proof of type _ inversion lemma is then

destruct v; ae.
It should be noted that this tactic does not produce any proof term, so we must trust both the way it reconstructs a Why3 task from the current Coq goal, and theorem provers called after, the same way we trust these theorem provers when called directly from Why3.

The first significant result states that a well-typed term of type t, is evaluated to a value of type t, too. This lemma has to be proved by structural induction on t. Once again, instead of making proof entirely in Coq, we first apply Why3 transformation for structural induction, which produces 4 new sub-goals, corresponding to a standard induction schema. Among these sub-goals, 3 are automatically proved (See Figure 8), the only one that we have to prove with Coq concerns binary operations which is done in 7 lines: In total, this proof of preservation of typing by reduction is only three lines long in Coq! Notice that 5 of the 7 seconds (see Figure 8) are lost when waiting for Alt-Ergo failing on the 8th sub-goal.

Substitutions, fresh variables

For defining our weakest precondition calculus, we will need to substitute variables. We dedicate a specific theory to the operation of substitution and its properties, that are particularly expressed under hypothesis of freshness of variables. Naming and freshness problems are classically difficult points to handle when we formalize languages with binders, as shown by the famous POPLmark challenge [START_REF] Aydemir | Mechanized metatheory for the masses: The POPLmark challenge[END_REF]. We propose here a basic approach, that we believe is easy to deal with, and sufficient for our study. First, we define the substitution operation, where we do not care about the problem of variable capture. The only one operation that we are interested in is to substitute a mutable program variable by a logical variable.

( ** substitution of a mutable variable [x] by a logic variable [v] warning: proper behavior only guaranteed if [v]

→ Fand (msubst f1 x v) (msubst f2 x v) | Fnot f → Fnot (msubst f x v) | Fimplies f1 f2 → Fimplies (msubst f1 x v) (msubst f2 x v) | Flet y t f → Flet y (msubst _ term t x v) (msubst f x v) | Fforall y ty f → Fforall y ty (msubst f x v) end RR n°8185
Then, we define the notion of fresh variable, that is, does not appear in a term or a formula. These are natural recursive definitions.

( ** [fresh _ in _ term id t] is true when [id] does not occur in [t] * ) predicate fresh _ in _ term (id:ident) (t:term) = match t with

| Tvalue _ | Tderef _ → true | Tvar i → id = i | Tbin t1 _ t2 → fresh _ in _ term id t1 ∧ fresh _ in _ term id t2 end predicate fresh _ in _ fmla (id:ident) (f:fmla) = match f with | Fterm e → fresh _ in _ term id e | Fand f1 f2 | Fimplies f1 f2 → fresh _ in _ fmla id f1 ∧ fresh _ in _ fmla id f2 | Fnot f → fresh _ in _ fmla id f | Flet y t f → id = y ∧ fresh _ in _ term id t ∧ fresh _ in _ fmla id f | Fforall y ty f → id = y ∧ fresh _ in _ fmla id f end
For reasoning about substitutions in the remaining of our formalization, we need some general lemmas we present now. Unsurprisingly, these lemmas have freshness hypothesis on concerned variables.

The two first lemmas relate the evaluation of a term (resp. a formula) on which we apply a substitution, with the evaluation of this term (resp. formula) in a modified state. In other words we have the identity These lemmas are proved by structural induction on the term (resp. formula), thanks to the Why3 transformation.

t[x ← v] Σ,Π = t Σ[x←Π(v)],
All sub-goals are automatically discharged (see Figure 10), except two which must be proved in Coq (case Flet and Fforall). These proofs are simple:

destruct f; auto. simpl; ae.

The next lemmas allow us, when evaluating a term or a formula, to swap two consecutive identifiers in the stack

(if they are different), that is t Σ,Π1•(id1,v1)•(id2,v2)•Π2 = t Σ,Π1•(id2,v2)•(id1,v1)•Π2
. The first of these lemmas, for terms, is written as6 lemma eval _ swap _ term: forall t:term, sigma:env, pi l:stack, id1 id2:ident, v1 v2:value. id1 = id2 → eval _ term sigma (l++(Cons (id1,v1) (Cons (id2,v2) pi))) t = eval _ term sigma (l++(Cons (id2,v2) (Cons (id1,v1) pi))) t

and is, again, proved by structural induction (See Figure 11). There is only one sub-goal which is not automatically discharged, the one corresponding to the case Tvar. This case is actually subtle, because it requires another induction, on the list l. In Coq, we use induction l tactic, and we finish the 2 sub-goals with why3 tactic.

Proof obligations

Alt-Ergo (0.93. We also pose the following instance of eval _ swap _ gen for the case l=Nil.

lemma eval _ swap: forall f:fmla, sigma:env, pi:stack, id1 id2:ident, v1 v2:value.

id1 = id2 → (eval _ fmla sigma (Cons (id1,v1) (Cons (id2,v2) pi)) f ↔ eval _ fmla sigma (Cons (id2,v2) (Cons (id1,v1) pi)) f)
While it is sufficient to instantiate the preceding lemma, this lemma is not proved by any theorem prover (see Figure 13). It has to be done in Coq intros f sigma pi id1 id2 v1 v2 h1. apply eval _ swap _ gen with (l:=Nil); auto.

Finally, two last lemmas allow to pop from the stack an identifier that does not appear in the evaluated term (resp. formula), that is t Σ,(id,v)•Π = t Σ,Π if id is fresh. lemma eval _ term _ change _ free : forall t:term, sigma:env, pi:stack, id:ident, v:value. fresh _ in _ term id t → eval _ term sigma (Cons (id,v) pi) t = eval _ term sigma pi t lemma eval _ change _ free : forall f:fmla, sigma:env, pi:stack, id:ident, v:value. fresh _ in _ fmla id f → (eval _ fmla sigma (Cons (id,v) pi) f ↔ eval _ fmla sigma pi f) Once again, the proof is done by structural induction (See Figure 14). The only sub-case which is not automatically proved is the case of Fforall in the direct way of implication, that is proved with: destruct f; auto. simpl; intros H sigma pi id v (h1 & h2). destruct d; intros; rewrite <-(H _ _ id v); ae.

Hoare Logic

This section is in fact not needed for the remaining of the case study. We just illustrate how one can formalize the standard Hoare logic rules. The notion of partial correctness of a Hoare triple, usually denoted {p}s{q}, is defined as follows.

predicate valid _ triple (p:fmla) (s:stmt) (q:fmla) = forall sigma:env, pi:stack. eval _ fmla sigma pi p → forall sigma':env, pi':stack, n:int. many _ steps sigma pi s sigma' pi' Sskip n → eval _ fmla sigma' pi' q

RR n°8185

Proof obligations

Alt-Ergo (0.93.

Figure 14: Substitution and freshness First, the consequence rule below can be proved automatically.

lemma consequence _ rule: forall p p' q q':fmla, s:stmt.

valid _ fmla (Fimplies p' p) → valid _ triple p s q → valid _ fmla (Fimplies q q') → valid _ triple p' s q'

RR n°8185

The rule for sequence is also proved automatically, but only after posing the standard "sequence lemma": lemma many _ steps _ seq: forall sigma1 sigma3:env, pi1 pi3:stack, s1 s2:stmt, n:int. many _ steps sigma1 pi1 (Sseq s1 s2) sigma3 pi3 Sskip n → exists sigma2:env, pi2:stack, n1 n2:int. many _ steps sigma1 pi1 s1 sigma2 pi2 Sskip n1 ∧ many _ steps sigma2 pi2 s2 sigma3 pi3 Sskip n2 ∧ n = 1 + n1 + n2 lemma seq _ rule: forall p q r:fmla, s1 s2:stmt.

valid _ triple p s1 r ∧ valid _ triple r s2 q → valid _ triple p (Sseq s1 s2) q

The proof of that lemma, even on paper, is not completely trivial and deserves a few lines, with an induction on the number n of steps. The proof is thus done in Coq, and is not surprising: In all these five cases, Coq must be used, and are done quite smoothly by inversion of predicates many _ steps and one _ step, and finishing the sub-goals with the Why3 tactic.

Theorem skip _ rule : forall (q:fmla), (valid _ triple q Sskip q). intros q sigma pi. intros H sigma' pi' n H1. inversion H1; subst; clear H1; auto. inversion H0. Qed.

Theorem assign _ rule : forall (p:fmla) (x:mident) (id:ident) (t:term), (fresh _ in _ fmla id p) → (valid _ triple (Flet id t (msubst p x id))

(Sassign x t) p). intros p x id t h1. red; intros. inversion H0; subst; clear H0. inversion H1; subst; clear H1. inversion H2; subst; clear H2. ae. ae. Qed.

Theorem if _ rule : forall (t:term) (p:fmla) (q:fmla) (s1:stmt) (s2:stmt), ((valid _ triple (Fand p (Fterm t)) s1 q) ∧ (valid _ triple (Fand p (Fnot (Fterm t))) s2 q)) → (valid _ triple p (Sif t s1 s2) q). unfold valid _ triple. intros t p q s1 s2 (h1,h2

→ q | Sassert f → Fand f (Fimplies f q)
( ** asymmetric and * ) | Sseq s1 s2 → wp s1 (wp s2 q) | Sassign x t → let id = fresh _ from q in Flet id t (msubst q x id)

| Sif t s1 s2 → ( * (t → WP(s1,Q)) ∧ (not t → WP(s2,Q)) * ) Fand (Fimplies (Fterm t) (wp s1 q)) (Fimplies (Fnot (Fterm t)) (wp s2 q)) | Swhile cond inv body → ( * inv ∧ forall effects, (cond ∧ inv → WP(body,inv)) ∧ (not cond ∧ inv → Q) * )
Fand inv (abstract _ effects body (Fand (Fimplies (Fand (Fterm cond) inv) (wp body inv)) (Fimplies (Fand (Fnot (Fterm cond)) inv) q))) end Our definition above follows the classical schema for such a calculus. For the case of Sassert f, we use f ∧ (f → Q) instead of f ∧ Q as it is done in practice 7 . For an assignment Sassign x t, we substitute the mutable variable x by a fresh logic variable id, computed thanks to a auxiliary function fresh _ from, then we bind id to the value assigned t. The case of the loop classically introduces a quantification on the side effects of the body of that loop. This quantification is realized by calling another auxiliary function abstract _ effects.

In order to deal with difficulties separately, we do not detail implementations of those auxiliary functions, we only axiomatize them. The first function is axiomatized by: function fresh _ from (f:fmla) : ident axiom fresh _ from _ fmla: forall f:fmla. fresh _ in _ fmla (fresh _ from f) f

In other words, for any formula f , fresh _ from f returns an arbitrary identifier that is fresh in f ; we don't care about how it is computed.

The function that abstract the side effects has type function abstract _ effects (s:stmt) (f:fmla) : fmla and its expected behavior is, for given statement s and formula P , to return the formula ∀w, f where w = w 1 , . . . , w k is the set of assigned variables in s. Our axiomatization of this behavior is made of four hypotheses detailed below. The first hypothesis formalizes the specialization property: if Σ, Π |= ∀w, f then Σ, Π |= f . axiom abstract _ effects _ specialize : forall sigma:env, pi:stack, s:stmt, f:fmla.

eval _ fmla sigma pi (abstract _ effects s f) → eval _ fmla sigma pi f

The second hypothesis formalizes the distributivity of quantification over conjunction : if Σ, Π |= (∀w, P )∧(∀w, Q) then Σ, Π |= ∀w, P ∧ Q.

axiom abstract _ effects _ distrib _ conj : forall s:stmt, p q:fmla, sigma:env, pi:stack.

eval _ fmla sigma pi (abstract _ effects s p) ∧ eval _ fmla sigma pi (abstract _ effects s q) → eval _ fmla sigma pi (abstract _ effects s (Fand p q))

The third hypothesis is a bit more subtle: it expresses a property of the quantification on the two parts of an implication.

axiom abstract _ effects _ monotonic : forall s:stmt, p q:fmla. valid _ fmla (Fimplies p q) → forall sigma:env, pi:stack. eval _ fmla sigma pi (abstract _ effects s p) → eval _ fmla sigma pi (abstract _ effects s q)

In other words, if |= P → Q then |= (∀w, P ) → (∀w, Q). This property can be seen as a consequence of two other properties: first, if |= f then |= ∀w, f ; second, if ∀w, (P → Q) then (∀w, P ) → (∀w, Q). It is important to notice that this hypothesis talks about validity in all states, because for fixed Σ and Π, Σ, Π |= P → Q does not imply Σ, Π |= (∀w, P ) → (∀w, Q). A counter-example is as follows: assuming Σ such that Σ(x) = 42, then true → x = 42 is valid, but not (∀x, true) → (∀x, x = 42).

Until now, the hypotheses on abstract _ effects express abstractly that it behaves as an universal quantification, without saying precisely on what. The fourth and last hypothesis we pose now expresses on what variables we quantify: if w = w 1 , . . . , w k is the set of variables assigned by s, then ∀w, P is a formula that is invariant by the WP calculus through s. axiom abstract _ effects _ writes : forall sigma:env, pi:stack, s:stmt, q:fmla. eval _ fmla sigma pi (abstract _ effects s q) → eval _ fmla sigma pi (wp s (abstract _ effects s q))

In fact, this abstract way of specifying the behavior abstract _ effects is precisely the important property that one tries to obtain when introducing the universal quantification into the weakest precondition of a loop: the formula cond ∧ inv → WP(body, inv) must be quantified, so that it becomes true (because it is invariant) at each iteration of the loop.

Basic Properties of the Calculus

In the literature, there exist advanced theoretical studies on the family of predicate transformers in general, that includes the weakest precondition calculus [START_REF] Back | Refinement calculus -a systematic introduction[END_REF]. Several properties are stated, often without showing why they are useful. In our case study, we are going to state and prove two properties that are needed for the main proof of soundness.

The first of these two properties is classically named monotonicity of the calculus. It is stated mathematically under the form: for any statement s and any formulas P and Q, if |= P → Q then |= WP(s, P ) → WP(s, Q). In Why3 it becomes:

Proof obligations

Alt-Ergo (0.93. valid _ fmla (Fimplies p q) → valid _ fmla (Fimplies (wp s p) (wp s q))

Notice that, as before, quantification on all states is essential: Σ, Π |= P → Q does not necessarily imply Σ, Π |= WP(s, P ) → WP(s, Q). There is a similar counter-example: if Σ(x) = 42 then (true → x = 42) but WP(x := 7, true) = true does not imply WP(x := 7, x = 42) = (7 = 42).

The proof of the monotonicity property is done by induction on the structure of s. In Why3, after application of the transformation for structural induction, it gives 6 sub-goals, as shown on Figure 16.

Only the cases of the assignment, the conditional and the loop are not proved automatically. The proofs for the case of Sassigns and Sif can both be done using the 2 same lines: destruct s; auto. unfold valid _ fmla; simpl; ae.

it should be noted that that last call to the why3 tactic makes use of lemmas for the substitution and fresh variables. The case of Swhile est slightly more involved, one needs to apply the hypothesis (3): destruct s; auto. unfold valid _ fmla; simpl. intros H1 p q H2; intuition. apply abstract _ effects _ monotonic with (2:=H3). unfold valid _ fmla; simpl. intuition.

The second property that we need is the distributivity over the conjunction: if Σ, Π |= WP(s, P ) and Σ, Π |= WP(s, Q) then Σ, Π |= WP(s, P ∧ Q). lemma distrib _ conj: forall s:stmt, sigma:env, pi:stack, p q:fmla. eval _ fmla sigma pi (wp s p) ∧ eval _ fmla sigma pi (wp s q) → eval _ fmla sigma pi (wp s (Fand p q))

The reciprocal is also true, but not needed. The proof proceeds again by structural induction on s, as shown on Figure 17. The cases of assignment, sequence, and loop are not proved automatically, the proofs are terminated in Coq.

The proof for the case of assignment is again easy: The proof in the case of the sequence is more complex. We need to use the preceding property of monotonicity. It indeed corresponds to an "intelligent" part of the proof, where we state as an intermediate lemma that

|= WP(s 1 , WP(s 2 , P ) ∧ WP(s 2 , Q)) → WP(s 1 , WP(s 2 , P ∧ Q)).
destruct s; auto. simpl; intros H1 H2 sigma pi p q (H3 & H4). assert (H: valid _ fmla (Fimplies (Fand (wp s2 p) (wp s2 q)) (wp s2 (Fand p q)))). unfold valid _ fmla ; simpl; ae. generalize (monotonicity s1 _ _ H). unfold valid _ fmla; simpl; ae.

In the case of the loop, we use hypotheses (3) and (4). destruct s; auto. simpl. intros H sigma pi p q ((h1 & h2) & (h3 & h4)). split; auto. apply abstract _ effects _ monotonic with (p:= (Fand (Fand (Fimplies (Fand (Fterm t) f) (wp s f))

(Fimplies (Fand (Fnot (Fterm t)) f) p)) (Fand (Fimplies (Fand (Fterm t) f) (wp s f))

(Fimplies (Fand (Fnot (Fterm t)) f) q)))). unfold valid _ fmla; simpl. intuition. apply abstract _ effects _ distrib _ conj; auto.

Soundness of the Weakest Precondition Calculus

We follow the classical approach for proving type soundness. It amounts to state two lemmas: on the one hand, the preservation of WP by reduction ; on the other hand a progress property, that is the validity of WP(s, Q) entails the reductibility of s (if s = Sskip). lemma wp _ preserved _ by _ reduction: forall sigma sigma':env, pi pi':stack, s s':stmt. one _ step sigma pi s sigma' pi' s' → forall q:fmla. eval _ fmla sigma pi (wp s q) → eval _ fmla sigma' pi' (wp s' q) RR n°8185

Proof obligations

Coq (8.3pl4) lemma wp _ preserved _ by _ reduction 6.78 (5s) (5s) (5s) (5s) 0.75 (5s) (5s)

Figure 19: Progress lemma This first lemma is proved by induction on hypothesis (one _ step sigma pi s sigma' pi' s'). This must be done in Coq (Figure 18). After using the induction tactic in Coq, there are 8 sub-goals, among which 6 can be proved simply by the simpl tactic followed by the why3 tactic. The two remaining cases are for the reduction of loops. To prove them, we start by using the hypothesis (1) to make explicit the fact that the formula (cond ∧ inv → WP(body, inv)) ∧ (¬cond ∧ inv → Q) is true in the current state, then we conclude using the why3 tactic. One should notice that these calls to why3 tactic use the hypothesis (4) and the property distrib _ conj. The progress lemma is stated as follows.

lemma progress: forall s:stmt, sigma:env, pi:stack, sigmat: type _ env, pit: type _ stack, q:fmla.

compatible _ env sigma sigmat pi pit ∧ type _ stmt sigmat pit s ∧ eval _ fmla sigma pi (wp s q) ∧ s = Sskip → reducible sigma pi s

It may seem strange that the lemma above depends on an arbitrary formula Q, that does not play any role in the conclusion. Thus, it may seem enough to state this lemma for Q = true. But then it would not be general enough to be proved by induction. For example, to prove the progress of a sequence s 1 ; s 2 knowing that WP(s 1 ; s 2 , true) holds, we need to prove the progress of s 1 knowing that WP(s 1 , WP(s 2 , true)) holds, and the formula WP(s 2 , true) is a priori arbitrary. The progress lemma is proved by structural induction on s. Only the case of Sskip is then proved automatically (Figure 19), the reason being that the conclusion reducible is existentially quantified. In each case we terminate the proof in Coq, in a few lines, with the exists tactic. For the case of the loop, we use once more the property Our main result is now: for any program s, any state Σ, Π and any formula Q, if Σ, Π |= WP(s, Q) then either s executes infinitely, or it reduces to Sskip, and in the second case Q is true in the final state. Since it is not handy to formulate that an execution is infinite, we state that slightly differently: if s executes into some finite number of steps into s such that s does not reduce anymore, then s = Sskip. All this must be stated only for well-typed programs, thus the Why3 formulation of our main result is: lemma wp _ soundness: forall n :int, sigma sigma':env, pi pi':stack, s s':stmt, sigmat: type _ env, pit: type _ stack, q:fmla. compatible _ env sigma sigmat pi pit ∧ type _ stmt sigmat pit s ∧ many _ steps sigma pi s sigma' pi' s' n ∧ not (reducible sigma' pi' s') ∧ eval _ fmla sigma pi (wp s q) → s' = Sskip ∧ eval _ fmla sigma' pi' q This is proved by induction on n, what we do calling Coq (Figure 20). For n = 0, the progress lemma ensures that s can only be Sskip. For the case n > 0, we know that s reduces in one step into some s 0 , the lemma of preservation of WP by reduction allows us to apply the induction hypothesis on s 0 .

intros n sigma sigma' pi pi' s s' sigmat pit q. intros Hcomp Htype (h1,(h2,h3)). generalize (steps _ non _ neg _ _ _ _ _ _ _ h1). intros Hnpos. generalize sigma sigma' pi pi' s s' sigmat pit q Hcomp Htype h1 h2 h3. clear sigma sigma' pi pi' s s' sigmat pit q Hcomp Htype h1 h2

General remarks on proofs

Let's put aside the section on Hoare logic which was presented for illustration but not needed for the main results. Overall, the proofs that we needed to do in Coq amount to 142 lines of tactics. This number should be compared with the number of lines of specification that we wrote: 390, that is more than twice the number of lines of proofs. This is a really good level of automation, that invalidates the usual belief that a proof of complex behavior of a program requires significantly more lines of proofs than lines of code.

The source of this case study, together with the Why3 proof session and the Coq proofs are available on the web page http://toccata.lri.fr/gallery/WP _ revisited.en.html. A ZIP archive is provided, that contains everything needed to replay the proofs using Why3 0.80 [START_REF] Bobot | The Why3 platform, version 0.80[END_REF].

Conclusions

We have presented a technique for proving the soundness of a weakest preconditions calculus, that imitates the classical technique of proof of type soundness of a programming language, by showing both that the WP is preserved by reduction and that the validity of the WP entails the reducibility. Since the conformance of the annotations of a program is by definition the absence of blocking (blocking semantics), we obtain a verification method, guaranteed correct, even on non-terminating programs. The definition via blocking semantics is proposed by Herms et al. [START_REF] Herms | A certified multi-prover verification condition generator[END_REF] in the case of a great-step semantics, soundness being proved by co-induction. We use a small-step blocking semantics : our approach is thus original as far as we know. Indeed, a similar approach was used by Conchon and Filliâtre [START_REF] Conchon | Semi-Persistent Data Structures[END_REF] in a different context: proof of soundness of a decision procedure for the safety of use of semi-persistent data structures.

The proof obtained by our approach is not difficult, in fact we pretend this is a natural approach, simple to understand and moreover is done with highly automated proof. We plan to use this method for teaching the basics of deductive program verification 8 .

Along this case study, we identified several interesting features of the Why3 environment: structuring in theories, algebraic data types, inductive predicate, direct calls to automated provers, transformation for proving by structural induction, call to Coq to terminate difficult proofs, use of the why3 tactic to terminate the Coq proofs. At the end, the number of lines of manual proof is small: 142 lines of Coq proofs for 390 lines of program code, that is a very satisfactory degree of automation. Nevertheless, one should be aware that during the development of this study, we had to do, temporarily, more Coq proofs, in order to identify the appropriate lemmas needed to finally get good automation. The why3 tactic is by the way very useful during such a phase of development of proofs.

We also identified a few weaknesses of Why3 that could make interesting motivations for extensions. The automated provers appear to be weak as soon as a hypothesis involves an inductive predicate: a transformation to make an inversion of such a hypothesis, so as to reason on it by induction, could be very handy and useful. For example, the lemma steps _ non _ neg could be proved without the need of Coq. Similarly, a transformation for proving a lemma by induction on integers would be useful.

We noticed a general weakness of SMT provers: often they can't prove a goal that can be proved in a few lines of Coq. For example, these provers are not able to do a simple thing like eexists; eauto of Coq. We also tested some provers specialized for first-order reasoning (of the kind « TPTP » : Vampire [START_REF] Riazanov | Vampire[END_REF], Spass [START_REF] Weidenbach | Spass version 3.5[END_REF], Eprover [START_REF] Schulz | System description: E 0.81[END_REF]) without success. We also noticed that the automated provers do not handle well the goals that can be proved very easily with the simpl tactic of Coq followed by a call to the why3 tactic. This emphasizes a weakness in the support for computations within proofs by such provers.

Another kind of missing feature is the ability to pose statements about the functions of the Why3 programming language. In fact, in this case study we only defined pure functions. If a function like wp was written using the programming language (that allow side-effects), then it would not be possible to state properties of its behavior like we did, but only by giving a post-condition, and in particular it would not be possible to state a property like wp _ soundness.

The case study itself deserves to be extended. First, the functions fresh _ from and abstract _ effects should be realized. A significant extension would to handle a language with several sub-programs, like it was done by Herms et al. [START_REF] Herms | A certified multi-prover verification condition generator[END_REF] and was indeed a reason to use the blocking semantics approach.

In a longer term, one could seriously think about using this approach to develop code correct by construction. In fact, an extension of Why3 currently in progress is the ability to extract OCaml code. With respect to an approach fully in Coq, we would gain a lot of automation.

1Figure 1 :

 1 Figure 1: Structure of the development

  type fmla = ( ** Formulas * ) | Fterm term ( ** terms: atomic formulas * ) | Fand fmla fmla ( ** conjunction * ) | Fnot fmla ( ** negation * ) | Fimplies fmla fmla ( ** implication * ) | Flet ident term fmla ( ** local binding: let id = term in fmla * ) | Fforall ident datatype fmla ( ** universal quantification: forall id : ty, fmla * ) Finally, the syntax of statements of our imperative language follows. type stmt = ( ** Statements * ) | Sskip ( ** no-op statement * ) | Sassign mident term ( ** assignment id := term * ) | Sseq stmt stmt ( ** sequence * ) | Sif term stmt stmt ( ** conditional statement * ) | Sassert fmla ( ** assertion statement * ) | Swhile term fmla stmt ( ** while loop with condition, invariant and body * )

Figure 2 :Figure 3 :

 23 Figure 2: One-step reduction

  inductive type _ operator (op:operator) (ty1 ty2 ty: datatype) = | Type _ plus : type _ operator Oplus TYint TYint TYint | Type _ minus : type _ operator Ominus TYint TYint TYint | Type _ mult : type _ operator Omult TYint TYint TYint | Type _ le : type _ operator Ole TYint TYint TYbool We define typing environments similarly to evaluation ones. type type _ stack = list (ident, datatype) ( ** map local immutable variables to their type * ) function get _ vartype (i:ident) (pi:type _ stack) : datatype = match pi with | Nil → TYunit | Cons (x,ty) r → if x=i then ty else get _ vartype i r end type type _ env = IdMap.map mident datatype ( ** map global mutable variables to their type * ) function get _ reftype (i:mident) (e:type _ env) : datatype = IdMap.get e i Typing judgments for terms, formulas and statements are, then, naturally defined by new inductive predicates, given on Figures 5, 6 and 7 inductive type _ term type _ env type _ stack term datatype = | Type _ value : forall sigma: type _ env, pi:type _ stack, v:value. type _ term sigma pi (Tvalue v) (type _ value v) | Type _ var : forall sigma: type _ env, pi:type _ stack, v: ident, ty:datatype. (get _ vartype v pi = ty) → type _ term sigma pi (Tvar v) ty | Type _ deref : forall sigma: type _ env, pi:type _ stack, v: mident, ty:datatype. (get _ reftype v sigma = ty) → type _ term sigma pi (Tderef v) ty | Type _ bin : forall sigma: type _ env, pi:type _ stack, t1 t2 : term, op:operator, ty1 ty2 ty:datatype.type _ term sigma pi t1 ty1 → type _ term sigma pi t2 ty2 → type _ operator op ty1 ty2 ty → type _ term sigma pi (Tbin t1 op t2) ty

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Typing rules for terms

  is fresh * ) function msubst _ term (t:term) (x:mident) (v:ident) : term = match t with | Tvalue _ | Tvar _ → t | Tderef y → if x = y then Tvar v else t | Tbin t1 op t2 → Tbin (msubst _ term t1 x v) op (msubst _ term t2 x v) end function msubst (f:fmla) (x:mident) (v:ident) : fmla = match f with | Fterm e → Fterm (msubst _ term e x v) | Fand f1 f2

Π

  lemma eval _ msubst _ term: forall e:term, sigma:env, pi:stack, x:mident, v:ident. fresh _ in _ term v e → eval _ term sigma pi (msubst _ term e x v) = eval _ term (IdMap.set sigma x (get _ stack v pi)) pi e lemma eval _ msubst: forall f:fmla, sigma:env, pi:stack, x:mident, v:ident. fresh _ in _ fmla v f → (eval _ fmla sigma pi (msubst f x v) ↔ eval _ fmla (IdMap.set sigma x (get _ stack v pi)) pi f)

Figure 18 :

 18 Figure 18: Preservation by reduction Proof obligations Alt-Ergo (0.93.1) Alt-Ergo (0.94) CVC3 (2.2) CVC3 (2.4.1) Coq (8.3pl4) Z3 (3.2) Z3 (4.2)

  Type env denotes environments for mutable variables. It is a map that associates a value to any mident. For this we use the Why3 standard library generic type map.

use export Syntax use map.Map as IdMap type env = IdMap.map mident value ( ** map global mutable variables to their value * ) function get _ env (i:mident) (sigma:env) : value = IdMap.get sigma i get _ env function provides a shortcut for accessing to the elements of an environment. Type stack denotes environments for logic variables: a stack of pairs (ident, value). We use the generic type list of Why3 standard library. use export list.List type stack = list (ident, value) ( ** map local immutable variables to their value * ) function get _ stack (i:ident) (pi:stack) : value = match pi with | Nil → Vvoid | Cons (x,v) r → if x=i then v else get _ stack i r end

  This lemma is proved in Coq, by induction on one _ step hypothesis. We obtain 8 sub-cases, 7 of them are automatically discharged by why3 tactic. The last sub-goal, the sequence, requires the hypothesis of well-typedness of the sequence before calling why3 tactic.

	destruct t; auto.
	simpl; intros.
	inversion H2; subst; clear H2.
	destruct H9 as (h1 & h2 & h3).
	generalize (type _ inversion (eval _ term sigma pi t1)).
	generalize (type _ inversion (eval _ term sigma pi t2)).
	destruct h3; ae.
	The second result is preservation of well-typing by reduction.
	lemma type _ preservation : forall s1 s2:stmt, sigma1 sigma2:env, pi1 pi2:stack,
	sigmat:type _ env, pit:type _ stack.
	type _ stmt sigmat pit s1 ∧ compatible _ env sigma1 sigmat pi1 pit ∧
	one _ step sigma1 pi1 s1 sigma2 pi2 s2 →
	type _ stmt sigmat pit s2 ∧ compatible _ env sigma2 sigmat pi2 pit
	intros s1 s2 sigma1 sigma2 pi1 pi2 sigmat pit (h1,(h2,h3)).
	induction h3; try ae.
	inversion h1; subst; clear h1; ae.

RR n°8185

http://why3.lri.fr/jfla-2012/ RR n°8185

Note that True and False are the constructors of Boolean type of Why3.RR n°8185

We remind that Why3 functions are always total.

The notion of inductive predicate in Why3 is similar to Coq or PVS ones: It is the smallest relation verifying the given clauses. Why3 verifies the usual monotonicity conditions for such definitions.RR n°8185

Indeed, we realized later in the development that CVC3 2.2 proves that lemma RR n°8185

Symbol ++ denotes concatenation of lists.RR n°8185

It forces f to appear as hypothesis in the other proof obligations.RR n°8185

rue Honoré d'Estienne d'Orves Bâtiment Alan Turing Campus de l'École Polytechnique 91120 Palaiseau Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments The authors thank Jean-Christophe Filliâtre and Paolo Herms for the concept of blocking semantics to define the validity of annotated code, Levs Gondelmans who implemented the transformation of structural induction of Why3, Jean-Christophe Filliâtre and Andrei Paskevich who implemented the why3 tactic, and also the other authors of Why3, Guillaume Melquiond and François Bobot, without whom this work would not have been achieved.

Cette étude a été partiellement financée par le projet U3CAT (ANR-08-SEGI-021, http://frama-c.com/u3cat/) de l'Agence Nationale de la Recherche et le projet Hi-lite (http://www.open-do.org/projects/hi-lite/) du pôle de compétitivité Systematic de la région Îlede-France.