Aparna Srivastav

Thomas Podgorski

Gwennou Coupier
email: gwennou.coupier@ujf-grenoble.fr

Efficiency of size-dependent particle separation by pinched flow fractionation

Keywords:

Pinched flow fractionation is shown to be an efficient and selective way to quickly separate particles by size in a very polydisperse semi-concentrated suspension. In an effort to optimize the method, we discuss the quantitative influence of the pinching intensity in the balance between the requirements of selectivity and minimal dilution.

Introduction

Sorting micro-particles of different sizes for analytic and preparative purposes has become a great challenge in the fields of chemical or biology research, or even industrial production. The development of microfluidic techniques has enabled control of flows at a scale similar to the size of the particles together with the use of small amounts of fluid, a target in the process of improving efficiency and reducing costs and dilution of samples.

As reviewed in [START_REF] Pamme | Continuous flow separations in microfluidic devices[END_REF]; Kersaudy-Kerhoas et al ( 2008)), several systems were recently developed in that purpose. Generally one expects a quick and precise separation of large quantities of particles, together with low cost and human intervention. In addition, the ability to handle highly concentrated suspensions can also be considered as a requisite. Among these systems, pinched flow fractionation (PFF), which was initially proposed in [START_REF] Yamada | Pinched flow fractionation: Continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF], has the advantage of Fig. 1 Scheme of the sorting device. The particles suspension is injected with flow rate Q 1 and is pinched by a particlefree fluid with flow rate Q 2 in a segment of width 2d. The suspension then enters a broad segment and splits into 30 outlets. In our experiment, 2d = 70µm and the thickness of the channels is 99µm. Drainage with flow rate Q 3 was also added, but its effect is not discussed here.

being a continuous process based only on hard-core interaction between particles and walls [START_REF] Luo | Irreversibility and pinching in deterministic particle separation[END_REF]) which is optimized by dedicated flow control. In particular, no use of an external field such as gravity or pressure waves is needed, so that no specific particle property is required.

The principle of PFF is extremely simple (see Figure 1): the suspension of (spherical) particles to be sorted is pushed against a wall by a particle-free pinching fluid. Due to their finite size, the centers of the particles are then located on the flow streamline one particle radius distant from the wall. In principle, collecting each of these streamlines leads to collecting each subpopulation of the sample. In practice, two geometrical improvements are necessary. As most of the pinching fluid will remain particle-free, it is convenient to limit its volume. The pinching is thus realized in a channel (the pinched segment) whose width 2d should be similar to the diameter 2R max of the largest particle to be sorted. Secondly, collecting subsamples with narrow size distributions requires that the lateral distance between the different populations is increased. This is achieved through a broadening of the pinched segment, allowing to add downstream collecting channels [START_REF] Andersen | A generalized theoretical model for "continuous particle separation in a microchannel having asymmetrically arranged multiple branches[END_REF]).

As at least the top half of the fluid will remain particle-free, it is convenient to drain it through a single outlet, in order to maximize the space available for the expansion of particle-charged streamlines [START_REF] Takagi | Continuous particle separation in a microchannel having asymmetrically arranged multiple branches[END_REF]; [START_REF] Vig | Separation enhancement in pinched flow fractionation[END_REF]; [START_REF] Maenaka | Continuous and size-dependent sorting of emulsion droplets using hydrodynamics in pinched microchannels[END_REF]). Alternatively, Sai et al. added microvalves to close some collecting outlets in order to better control the final destination of the particles [START_REF] Sai | Continuous separation of particles using a microfluidic device equipped with flow rate control valves[END_REF]). At low Reynolds number, which generally applies to such systems, the final destinations in the collecting set-up of the fluid streamlines created after pinching should theoretically not depend on the geometry of the pinched segment and its expansion. Some works however explored possible geometrical variations and their consequences on separation efficiency as compared to the expected one [START_REF] Yamada | Pinched flow fractionation: Continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF]; [START_REF] Zhang | Continuous flow separation of particles within an asymmetric microfluidic device[END_REF]; [START_REF] Maenaka | Continuous and size-dependent sorting of emulsion droplets using hydrodynamics in pinched microchannels[END_REF]). Anyhow, it is likely that the small variations observed are mainly due to channel imperfections, 3D effects or optical errors [START_REF] Jain | Particle dispersion and separation resolution of pinched flow fractionation[END_REF]).

In most of the above cited papers, the proposed setup is validated by injecting a very dilute mixture of a small number (between 2 and 4) of subpopulations of spherical particles with very different radii. Noting λ the ratio between the radii of a large particle and a smaller one, the following papers include such a validation: [START_REF] Yamada | Pinched flow fractionation: Continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF], λ = 2; [START_REF] Takagi | Continuous particle separation in a microchannel having asymmetrically arranged multiple branches[END_REF], λ ≥ 1.4; [START_REF] Zhang | Continuous flow separation of particles within an asymmetric microfluidic device[END_REF], λ = 2.5; [START_REF] Sai | Continuous separation of particles using a microfluidic device equipped with flow rate control valves[END_REF], λ ≥ 1.4; [START_REF] Jain | Particle dispersion and separation resolution of pinched flow fractionation[END_REF], λ = 1.5; [START_REF] Morijiri | Sedimentation pinched-flow fractionation for size-and density-based particle sorting in microchannels[END_REF], λ = 1.7. Vig and Kristensen considered seven subpopulations which are mixed and sorted, but their concentration is very low (0.05%) [START_REF] Vig | Separation enhancement in pinched flow fractionation[END_REF]). As a result, those studies do not fully validate all expected performances of a sorting device, namely the possibility to separate efficiently particles out of a continuum of sizes in a concentrated sample while avoiding excessive dilution in the process. As concentration may lead to increased interactions between particles and possibly decrease the selectivity of the system, the influence of this parameter needs to be explored. Maenaka et al. have considered an emulsion with continuous droplet radius distribution between 0 and 30 µm, with a concentration of 5% [START_REF] Maenaka | Continuous and size-dependent sorting of emulsion droplets using hydrodynamics in pinched microchannels[END_REF]). The sample is injected at a flow rate Q 1 = 300µL/h, with a confining flow Q 2 = 9Q 1 (see Figure 1), and three collecting channels are used. The size distributions in these three channels are quite well separated, although overlappings of around 5 microns seem to occur. As the initial population is only split into three sub-populations, the monodispersity, that is, the existence of a narrow size distribution relatively to the mean size, remains weak in each of them .

In this paper, we wish to go further and optimize the flow rate ratio Q 1 /Q 2 between the suspension and the pinching fluid, in order to find the good balance between the requirements of selectivity and minimal dilution. In the meantime, we show that subsamples with good monodispersity can be quickly obtained.

Experimental set-up

As test particles, we use lipid vesicles obtained through the standard electroformation technique [START_REF] Angelova | Preparation of giant vesicles by external ac electric fields. kinetics and applications[END_REF]), which straightforwardly produces polydisperse samples with a continuous size distribution. In order to ensure spherical shapes and prevent lateral migration in the pinched segment due to viscous forces [START_REF] Coupier | Noninertial lateral migration of vesicles in bounded Poiseuille flow[END_REF]), a sucrose solution is encapsulated and vesicles are immersed in a slightly hypo-osmotic glucose solution. Being spherical, the vesicles can only deform if their membrane is stretched. The ability of the flow to stretch the membrane is given by the capillary number C a = ηU/κ s , where η is the fluid viscosity, U its typical velocity, and κ the membrane stretching modulus, of order 0.3 N/m for DOPC membranes [START_REF] Rawicz | Effect of chain length and unsaturation on elasticity of lipid bilayers[END_REF]). The maximum flow rate considered here will be 20 mL/h, in a channel of typical size 100 × 100 µm, therefore C a will be of order 10 -3 at maximum, and vesicle deformation can be neglected. The use of different fluids inside and outside the vesicles also allows particle visualization through a phase contrast microscope coupled to a fast camera. The sorting device is a standard PDMS microfluidic chip bounded to a glass plate (Figure 1). The pinched segment has width 2d = 70µm and thickness 99µm. From now on, we set d as the lengthscale of the problem. The sorted subsamples are observed in 30 collecting channels located at the end of the broadened segment. In this proof of concept device with no sample collection, they all converge to a unique outlet at atmospheric pressure. Additional drainage through a sucking with flow rate Q 3 = 0.9(Q 1 + Q 2 ) was added in order to increase the number of channels with particles. Note that the best location of this drainage is a complex issue as an infinity of geometrical variations are possible [START_REF] Yamada | Pinched flow fractionation: Continuous size separation of particles utilizing a laminar flow profile in a pinched microchannel[END_REF]; [START_REF] Takagi | Continuous particle separation in a microchannel having asymmetrically arranged multiple branches[END_REF]; [START_REF] Vig | Separation enhancement in pinched flow fractionation[END_REF]; [START_REF] Maenaka | Continuous and size-dependent sorting of emulsion droplets using hydrodynamics in pinched microchannels[END_REF]). We shall not explore this in this work.

The particles have radii lying between 0 and 0.97 (in d unit), and split into the 16 first channels (which indicates that our drainage is not optimal). In each outlet channel the mean radius < R > i of the particles and the standard deviation is measured. The flow rate Q 2 was kept constant to 10 mL/h, and Q 1 was varied between 0.2 and 10 mL/h. We wish to find the optimum ratio Q 1 /Q 2 but this parameter would depend on the chosen width 2d (larger d would require larger Q 2 for the same pinching efficiency), which depends itself on the maximum size of the particles to be sorted. A more appropriate parameter is therefore the width d 1 occupied in the pinched-segment by the fluid coming from channel 1. Roughly, d 1 should be of the order of the radius of the smallest particles to be sorted.

Finally, initial volume concentrations of 0.8 and 4.8 % are considered.

Results

Figures 2(a) and 2(b) show the mean sizes of the particles in the different channels for different values of d 1 between 0.16 (Q 1 = 0.2 mL/h) and 1 (Q 1 = 10 mL/h). For d 1 ≤ 0.51, good separation is achieved, with standard deviations of the order of the half distance between two neighboring mean values.

In addition to this separation, the obtained subsamples have a good monodispersity: the particles of radii ranging from 0.05 to 0.97 split into 16 subsamples where the mean radius increases quasi linearly and the standard deviation is roughly constant in the different channels and equal to 0.017 for d 1 = 0.16, 0.021 for d 1 = 0.36 and to a still reasonable 0.036 for d 1 = 0.51. Note that if particles of radii between 0 and 1 are perfectly separated into 16 subpopulations, the size difference between the smallest and the largest particle in each subpopulation is 0.063, which gives a standard deviation of around 0.016, so with d 1 = 0.16 we reach the best possible monodispersity with our choice of outlet channels.

The monodispersity quality can be estimated by the ratio between the standard deviation and the mean radius, and from channel 1 to channel 16 we find that it goes from 17% to 2% when d 1 = 0.16 and from 36% to 4% when d 1 = 0.51. The low quality for small particles is related to the choice of outlet channels with equal widths, but the good quality obtained for larger particles show that this is just a scale issue that would be solved by local refinements at the level of the outlets receiving the smallest particles.

As d 1 is increased even more, the standard deviations increase drastically while the mean radius varies less between the outlets, in particular for high radii, so separation is not achieved. Note that, for d 1 = 0.62, outlets 14 to 16 have narrow standard deviations, as well as outlet 16 for d 1 = 0.77: the pinching is still

.8% d 1 =0.51 (Q 1 /Q 2 =0.2) d 1 =0.62 (Q 1 /Q 2 =0.3) d 1 =0.77 (Q 1 /Q 2 =0.5) d 1 =1.00 (Q 1 /Q 2 =1.0)
<R> Outlet d 1 =0.51 (Q 1 /Q 2 =0.2), concentration 0.8% d 1 =0.36 (Q 1 /Q 2 =0.1), concentration 4.8%
(c) concentration 0.8% d 1 =0.16 (Q 1 /Q 2 =0.02) d 1 =0.36 (Q 1 /Q 2 =0.1) d 1 =0.51 (Q 1 /Q 2 =0.2)
Fig. 2 (color online) Size distributions in the outlet channels. (a) and (b): results for initial concentration of 0.8% and varying pinched width d 1 (that is, varying Q 1 /Q 2 ). (c): results for initial concentration of 4.8% compared with previous ones.

strong enough to prevent smaller particles from being at the same level as the large particles that enter these outlets. For intermediate channels (around channel 10) or for d 1 = 1, the population is roughly bidisperse. This can be seen in Figure 3(a): the narrow distribution of outlet 3 becomes wider in outlet 7 and 10 and two distinct populations appear in outlet 12 and 15. A possible explanation of this phenomenon can be reached through purely geometrical considerations if one takes into account the effects of walls in inlet 1, before pinching. In this area, the distribution of small particles extends to closer distances to the top wall than larger ones. If the pinching is not strong enough, a reminiscence of this distribution will be found in the pinched fluid, with some small particles higher than larger ones. On the other hand, the bottom wall in the pinched segment prevent the largest particles from respecting this initial order and they must stay in upper position, where non-pinched small particles lie.

Finally, it is remarkable that separation with d 1 around 0.5 is almost as good as with the more intuitive choice of d 1 = 0.16, which is of the same order as the smallest particles considered here (even smaller ones, that enter channel 1, were barely distinguishable). If all the small particles present in the initial suspension were present in the whole pinched area of width d 1 = 0.51, they should be present until channel 10, where < R >= 0.50, which is clearly not the case. One must therefore admit that collective effects take place, so that for instance small particles could be pushed against the wall by larger particles.

Finally, as shown in Figures 2(c) and 3(b), good separation and monodispersity are also achieved for a concentration of 4.8%: standard deviations are 0.027 with d 1 = 0.36, which is comparable with the more dilute case. The monodispersity quality reaches 3% for the larger particles.

Conclusion

We have shown that pinched flow fractionation is an efficient technique to separate a semi-concentrated polydisperse suspension of micrometric spherical particles into subsamples with tiny overlappings. In addition, monodisperse suspensions with a few % of variation in sizes can be obtained from that initial suspension where size variations reach 50%. Moreover, broader initial size distribution would probably not increase the final monodispersity quality in each subsample since the main issue is to separate small to medium sized particles. According to the desired monodispersity quality of the final subsamples the location and width of the collecting outlets can be optimized. Here, we considered outlets of equal width, which resulted in constant standard deviation between subsamples, thus low relative monodispersity quality in the small particles samples. Getting constant monodispersity quality requires to consider outlet of linearly increasing width. If one wishes to get concentrated samples at the outlet, a pinched suspension width 4 or 5 times larger than the radius of the smallest particles surprisingly appears to be a good compromise: confining more tightens a little bit the distributions, but dilutes more the samples, while a weaker confinement leads to bad sorting. These conclusions are valid for a semi-concentrated suspension (concentrations around 5%) and for quite high flow rates (some mL/h, that is one order of magnitude higher than in [START_REF] Maenaka | Continuous and size-dependent sorting of emulsion droplets using hydrodynamics in pinched microchannels[END_REF]). Note that for Q 1 + Q 2 = 11 mL/h, the Reynolds number is of order 100 in the pinched segment, so lateral drift of inertial origin could have perturbed the sorting [START_REF] Segre | Radial particle displacements in poiseuille flow of suspensions[END_REF]). However, inertial lift or viscous lift (that can occur in the case of deformable particles) often increase with particle size, therefore they should preserve the sorting effect and only shift the location of collecting outlets.

In addition, we have shown that bidisperse suspensions can be easily obtained through weaker pinching.

Fig. 3 (

 3 Fig.3(color online) Particle probability densities n i (R) in the outlets. n i (R)dR is the probability to get a vesicle of radius between R and R + dR when picking up one in channel i. (a): for some selected outlets for a concentration of 0.8% and d 1 = 0.77 (same data as in 2(b)). (b): in the 14 outlets analyzed for a concentration of 4.8% and d 1 = 0.36 (same data as in 2(c)).


Arnaud Gotlieb
email: arnaud@simula.no

Tristan Denmat
email: tristan.denmat@inria.fr

Nadjib Lazaar
email: lazaar@lirmm.fr

Constraint-based reachability

Iterative imperative programs can be considered as infinite-state systems computing over possibly unbounded domains. Studying reachability in these systems is challenging as it requires to deal with an infinite number of states with standard backward or forward exploration strategies. An approach that we call Constraint-based reachability, is proposed to address reachability problems by exploring program states using a constraint model of the whole program. The keypoint of the approach is to interpret imperative constructions such as conditionals, loops, array and memory manipulations with the fundamental notion of constraint over a computational domain. By combining constraint filtering and abstraction techniques, Constraint-based reachability is able to solve reachability problems which are usually outside the scope of backward or forward exploration strategies. This paper proposes an interpretation of classical filtering consistencies used in Constraint Programming as abstract domain computations, and shows how this approach can be used to produce a constraint solver that efficiently generates solutions for reachability problems that are unsolvable by other approaches.

Introduction

Modern automated program verification can be seen as the convergence of three distinct approaches, namely Software Testing, Model-Checking and Program Proving. Even if the general verification problems are often undecidable, investigations on these approaches have delivered the most efficient automated techniques to show that a given property is satisifed or not by all the reachable states of an infinite-state system.

Several authors have advocated the usage of constraints to represent an infinite set of states and the usage of constraint solvers to efficiently address reachability problems [START_REF] Richard | Testing abstract distributed programs and their implementations: A constraintbased approach[END_REF][START_REF] Delzanno | Constraint-based deductive model checking[END_REF][START_REF] Flanagan | Automatic software model checking via constraint logic[END_REF][START_REF] Boonstoppel | RWset: Attacking path explosion in constraintbased test generation[END_REF]. In automated program verification problems, the goal is to find a state of the program which violates a given safety property, i.e., an unsafe state. Two distinct strategies have been investigated to explore programs with constraints, namely the forward analysis and the backward analysis strategies. In forward analysis, a set of reachable states is explored by computing the transition from the initial states of a program to the next states in forward way. If an unsafe state is detected to belong to the set of reachable states during this exploration then a property violation is reported. In backward analysis, states are computed from an hypotetical unsafe state in a backward way with the hope to discover that one of those is actually an initial state. One advantage of backward analysis over forward analysis is its usage of the targeted unsafe state to refine the state search space. However, both strategies are quite powerful and have been implemented into several software model checkers based on constraint solving [START_REF] Henzinger | Software verification with Blast[END_REF][START_REF] Flanagan | Automatic software model checking via constraint logic[END_REF] and automated test case generators [START_REF] Williams | PathCrawler: Automatic Generation of Path Tests by Combining Static and Dynamic Analysis[END_REF][START_REF] Godefroid | DART: directed automated random testing[END_REF][START_REF] Fraser | Test-Case Generation and Coverage Analysis for Nondeterministic Systems Using Model-Checkers[END_REF][START_REF] Charreteur | Constraint-Based Test Input Generation for Java Bytecode[END_REF][START_REF] Bardin | OSMOSE: Automatic Structural Testing of Executables[END_REF].

In this paper, we present an integrated constraint-based strategy that can benefit from the strengths of both forward and backward analysis. The keypoint of the approach, that we have called Constraint-Based Reachability (CBR), is to interpret imperative constructions such as conditionals, loops, array and memory manipulations with the fundamental notion of constraint over a computational domain. By combining constraint filtering and abstraction techniques, CBR is able to solve reachability problems which are usually outside the scope of backward or forward exploration strategies. A main difference is that CBR does not sequentially explore the execution paths of the program ; the exploration is driven by the constraint solver which picks-up the constraint to explore depending on the priorities that are attached to them. It is worth noticing that applying CBR to program exploration results in a semi-correct procedure only, meaning that there is no termination guarantee. CBR has been mainly applied in automatic test data generation for iterative programs [START_REF] Gotlieb | Automatic Test Data Generation Using Constraint Solving Techniques[END_REF][START_REF] Gotlieb | A CLP Framework for Computing Structural Test Data[END_REF], programs that manipulate pointers towards named locations of the memory [START_REF] Gotlieb | Constraint-based test data generation in the presence of stack-directed pointers[END_REF][START_REF] Gotlieb | Goal-oriented test data generation for pointer programs[END_REF], programs on dynamic data structures and anonymous locations [START_REF] Charreteur | Modelling dynamic memory management in Constraint-Based Testing[END_REF], programs containing floating-point computations [START_REF] Botella | Symbolic execution of floating-point computations[END_REF]. A major improvement of the approach was brought by the usage of Abstract Interpretation techniques to enrich the filtering capabilities of the constraints used to represent conditionals and loops [START_REF] Denmat | An Abstract Interpretation Based Combinator for Modeling While Loops in Constraint Programming[END_REF][START_REF] Denmat | Improving Constraint-Based Testing with Dynamic Linear Relaxations[END_REF]. This approach permitted us to build efficient test data generator tools for a subset of C [START_REF] Gotlieb | EUCLIDE: A Constraint-Based Testing platform for critical C programs[END_REF] and Java Bytecode [START_REF] Charreteur | Constraint-Based Test Input Generation for Java Bytecode[END_REF].

The first contribution of this paper is the interpretation of classical filtering consistencies notions in terms of abstract domain computations. Constraint filtering is the main approach behind the processing of constraints in a finite domains constraint solver. We show in general the existence of tight links between classical filtering techniques and abstract domain computations that were not pointed out elsewhere. We also give the definition of a new consistency filtering inspired from the Polyhedral abstract domain, as consequence of these links.

The second contribution is the description of a special constraint handling any iterative construction. The constraint w captures iterative reasoning in a constraint solver and as such, is able to deduce information which is outside the scope of any pure forward or backward abstract analyzer. Its filtering capabilities combines both constraint reasoning and abstract domain computations in order to propagate informations to the rest of the constraint system. In this paper, we focus on the theoretical foundations of the constraints, while giving examples of its usage for test case generation over iterative programs.

Outline of the paper. The rest of the paper is organized as follows. Sec.2 introduces the necessary background in Abstract Interpretation to understand the contributions of the paper. Sec.3 establishes the link between classical constraint filtering and abstract domain computations. Sec.4 describes the theoretical foundation of the w constraint for handling iterative constructions while Sec.5 concludes the paper.

Background

Abstract Interpretation (AI) is a theoretical framework introduced by Cousot and Cousot in [START_REF] Cousot | Abstract Interpretation : A unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] to manipulate abstractions of program states. An abstraction can be used to simplify program analysis problems otherwise not computable in realistic time, to manageable problems more easily solvable. Instead of working on the concrete semantics of a program1 , AI computes results over an abstract semantics allowing so to produce over-approximating properties of the concrete semantics. In the following we introduce the basic notions required to understand AI. Definition 1 (Partially ordered set (poset)) Let ⊑ be a partial order law, then the pair

(D, ⊑) is called a poset iff ∀x ∈ D, x ⊑ x (reflexive) ∀x, y ∈ D, x ⊑ y ∧ y ⊑ x =⇒ x = y (anti-symmetry) ∀x, y, z ∈ D, x ⊑ y ∧ y ⊑ z =⇒ x ⊑ z (transitive) Definition 2 (Complete lattice) A complete lattice is a 4-tuple (D, ⊑, , ) such that • (D, ⊑) is a poset • is a upper bound: ∀S ⊆ D, we have ∀x ∈ S , x ⊑ S ∀y ∈ D, (∀x ∈ S , x ⊑ y) =⇒ S ⊑ y • is a lower bound: ∀S ⊆ D, we have ∀x ∈ S , S ⊑ x ∀y ∈ D, (∀x ∈ S , y ⊑ x) =⇒ y ⊑ S
Complete lattices have a single smallest element ⊥ = D and a single greatest element ⊤ = D.

Program semantics can usually be expressed as the least fix point of a monotonic and continuous function.

A function f from a complete lattice (D, ⊑, ,

) to itself is monotonic iff ∀l 1 , l 2 ∈ D, l 1 ⊑ l 2 =⇒ f (l 1 ) ⊑ f (l 2 ). It is continuous iff ∀S ⊆ D, f ( S ) = s∈S ( f (s)) and f ( S ) = s∈S ( f (s)).
The following Theorem guarantees the existence of the fix points of a monotonic function.

Theorem 1 (Knaster-Tarski) In a complete lattice (D, ⊑, , ), for all monotonic functions f : D → D,

• the least fix point of f (i.e., l f p( f )) exists and l f p( f

) = {x | f (x) ⊑ x} • the greatest fix point of f (i.e., g f p( f )) exists and g f p( f ) = {x | f (x) ⊑ x}
In addition, when the functions are continuous, these fix points can be computed using an algorithm derived from the following theorem:

Theorem 2 (Kleene) In a complete lattice (D, ⊑, , ), for all monotonic and continuous functions f

: D → D, the least fix point of f is equal to { f n (⊥) | n ∈ N} and the greatest fix point of f is equal to { f n (⊤) | n ∈ N} As ⊥, f (⊥), . . . , f n (⊥), . . . is an increasing suite, we get { f n (⊥) | n ≤ k} = f k (⊥). Hence, lfp( f ) = lim k→+∞ f k (⊥) and gfp( f ) = lim k→+∞ f k (⊤).
For reaching the least fix point of a monotonic and continuous function in a complete lattice, it suffices to iterate f from ⊥ until a fix point is reached. Let (D, ⊑, , ) be a complete lattice called the concrete lattice and f a function that defines some concrete semantics over this lattice, let (D ♯ , ⊑ ♯ ) be a poset called the abstract poset, and f ♯ : D ♯ → D ♯ be a continuous function, then Abstract Interpretation aims at computing a fix point of f ♯ in order to over-approximate the computation performed by f . Depending on whether the abstract poset is a complete lattice or not, we have distinct theoretical results regarding the abstraction. Proofs of the following theorems can be found in [START_REF] Cousot | Automatic Discovery of Linear Restraints Among Variables of a Program[END_REF].

Galois connection

When the abstract poset is a complete lattice, the notion of Galois connection is available to link the abstract computations with the concrete lattice.

(D, ⊑, , ) γ ⇆ α (D ♯ , ⊑ ♯ , ♯ , ♯ )
Next definition establishes the correction property of an analysis.

Definition 4 (Sound approximation)

Let (D, ⊑, , ) γ ⇆ α (D ♯ , ⊑ ♯ , ♯ , ♯ ) be a Galois connection, then a function f ♯ : D ♯ → D ♯ is a sound approximation of f : D → D iff ∀y ∈ D ♯ , f • γ(y) ⊑ γ • f ♯ (y)
Consquently, we have the following notion: This theorem implies that any function greater than α • f • γ is a sound approximation of f and the following theorem characterizes the results of fixpoint computations:

Theorem 4 (Fixpoint computations with sound approximation) Let (D, ⊑, ,

) γ ⇆ α (D ♯ , ⊑ ♯ , ♯ , ♯ )
be a Galois connection, let f ♯ : D ♯ → D ♯ and f : D → D be two monotonic functions such that f ♯ is a sound approximation of f , then, we have:

lfp( f ) ⊑ γ(lfp( f ♯ )) and gfp( f ) ⊑ γ(gfp( f ♯ ))
Intuitively, this theorem gives a process to compute an over-approximation by Abstract Interpretation, as shown in Fig 1 . The left part shows the concrete lattice where the concrete computation of f is performed starting from initial state S 0 . The right part shows the abstract lattice that is used to over-approximate the computation. This computation is undertaken in three steps:

• initial state abstraction;

• fixpoint computation in the abstract lattice;

• result concretization.

Without Galois connection When the abstract lattice is not complete, there does not exist necessarily a best abstraction for all elements of the concrete lattice. The notion of Galois connection is no more available and the abstract lattice is just linked with the concrete lattice through a monotonic function γ : D ♯ → D. The definition of sound approximation needs to be adapted:

♯ : D ♯ → D ♯ is a sound approximation of f iff ∀x ∈ D ♯ , f • γ(x) ⊑ γ • f ♯ (x)
In such an (not complete) abstract lattice, nothing guarantees the existence of the least fix point: lfp( f ) is not necessarily approximated by lfp( f ♯ ). However, any fix point of f ♯ can be used:

Theorem 5 Let (D, ⊑, , ) be a complete lattice, and (D ♯ , ⊑ ♯ ) be a poset, let γ :

D ♯ → D, f : D → D and f ♯ : D ♯ → D ♯ be three monotonic functions then if f ♯ is a sound approximation of f , then we have: ∀x ∈ D ♯ , f ♯ (x) = x =⇒ lfp( f ) ⊑ γ(x)
Next theorem is useful to compute an over-approximation of gfp( f ) when the lattice is not complete:

Theorem 6 Let (D, ⊑, , ) be a complete lattice, let (D ♯ , ⊑ ♯ ) be a poset with a greatest element ⊤ and let γ : D ♯ → D, f : D → D and f ♯ : D ♯ → D ♯ be three monotonic functions, then if f ♯ is a sound approximation of f and a is an element of D ♯ such that there exists k such as a

= f ♯ k (⊤), then gfp( f ) ⊑ γ(a)
Consequently, when the abstract lattice is not complete, instead of abstracting the initial state, one selects an element of the abstract lattice that over-approximates the initial state. And, a fix point is computed in the abstract lattice from this element. The fix point is still an over-approximation of the concrete semantics.

Examples of abstract domains

In this section, we briefly describe two abstract domains: the Interval [START_REF] Cousot | Static Determination of dynamic properties of programs[END_REF] and the Polyhedral [START_REF] Cousot | Automatic Discovery of Linear Restraints Among Variables of a Program[END_REF] domains.

The Interval abstract domain

Interval analysis aims at approximating a set of values by an interval of possible values. If

I N = {[a, b] | a, b ∈ N ∪ {-∞,
α({(x 11 , . . . , x n1 ), (x 12 , . . . , x n2 ), . . .}) = ([m 1 , M 1 ], . . . , [m n , M n ]) with m i = min j (x i j ) if it exists -∞ else and M i = max j (x i j ) if it exists +∞ else
The concretization of an abstract state is obtained by computing the Cartesian product of the intervals. These functions define a Galois connection between the concrete domain and the abstract domain of intervals.

The approximation of transfert functions is realized by using their structure and classical results from Interval Analysis [START_REF] Moore | Interval Analysis[END_REF]. For example, functions [x ≤ y] and [x = y + z] are abstracted by the following (sound) approximations:

[x ≤ y] ♯ : ([a, b], [c, d]) → ([a, min(b, d)], [max(a, c), d])) and [x = y + z] ♯ : ([a, b], [c, d], [e, f ]) → ([c + e, d + f ] ∩ ♯ [a, b], [a -f , e -b] ∩ ♯ [c, d], [a -d, c -b] ∩ ♯ [e, f ]).

The Polyhedral abstract domain

In Polyhedral analyses, each concrete state is abstracted by a conjunction of linear constraints that defines a convex polyhedron. Indeed, a convex polyhedron is a region of an n-dimensional space that is bounded by a finite set of hyperplanes x ∈ IR n |ax ≥ c where a ∈ IR n and c ∈ IR. The abstract lattice equiped with inclusion, convex hull2 , and intersection of polyhedra is not a complete lattice as there is no upper bound to the convex union of all the convex polyhedra that can be written in a circle. Abstract functions can be defined to deal with polyhedra. For example:

[x ≥ y] ♯ ({z ≤ x + y}) = {z ≤ x + y, y ≤ x} (1) [x > y] ♯ ({x ≤ y}) = {0 = 1} (2) [x = y * z] ♯ ({1 ≤ y ≤ 10}) = {x ≤ z, x ≤ 10 * z} (3)
If the expression is a linear condition, then it is just added to the polyhedron (case 1). If the expression is contradictory with the current polyhedron, then it is reduced to 1 = 0 meaning that there is no abstract (and concrete) state in the approximation (case 2). If the expression is non-linear, then a linear approximation is derived when available and added to the polyhedron (case 3).

Filtering consistencies as abstract domain computations

As noticed by Apt [START_REF] Apt | The Essence of Constraint Propagation[END_REF], constraint propagation algorithms can be seen as instances of algorithms that deal with chaotic iteration. In this context, chaotic means fair application of propagators until saturation.

In this section, we elaborate on a bridge between two unrelated notions: filtering consistencies and abstract domains. In particular, we show that arc-and bound-consistency are instances of chaotic iterations over two distinct abstract domains. Classical AI notions of sound approximation and abstract domain computations, not used in [START_REF] Apt | The Essence of Constraint Propagation[END_REF], allows to show that filtering consistencies compute sound overapproximations of the solutions set of a constraint system. Thanks to the bridge, we also propose new filtering consistency algorithms based on the polyhedral abstract domain.

Notations

Let Z Z be the set of integers and V be a finite set of integer variables, where each variable x in V is associated with a finite domain D(x). The domain D is the Cartesian product of each variable domain: D(x 1 ) × . . . × D(x m ) and P(D) denotes the powerset of D. in f D x and sup D x denote respectivelly the inferior and the superior bounds of D(x) in D. A constraint c is a relation between variables of V . The language of (elementary) constraints is built over arithmetical operators {+, -, * , ...} and relational operators {<, ≤, >, ≥, =, =, ...} but any relation over a subset of V can be considered. Let vars(c) be the function that returns the variables of V appearing in a constraint c. A valuation σ is a mapping of variables to values, noted {x 1 → d 1 , ..., x n → d n }. CS denotes a constraint system CS, i.e., a finite set of constraints.

Exact filtering

Let {c 1 , .., c m } be a CS over {x 1 , .., x n } and let D = D(x 1 ) × .. × D(x n ), then the solution-set of CS is an element of P(D), noted sol(CS).

The exact filtering operator of a constraint c i is computed with the function f i : P(D) → P(D) which maps an element S ∈ P(D)

to f i (S) = {s | s ∈ S ∧ c i (s)}.
The exact filtering operator of c i removes all the tuples of D that violate c i . Hence, by using an iterating procedure, it permits to compute sol(CS):

if f C = f 1 • . . . • f m then sol(C) = gfp( f C )
. By noticing that f C is continuous (as each f i is continuous) and monotonic and thanks to Theorem 2 we get sol(CS)

= lim k→+∞ f C k (D). Example 1 Consider CS = {x = y, y = z, z = x} where x ∈ 1..2, y ∈ 1..2, z ∈ 1..2.
The exact filtering operator associated with x = y will remove the tuples (

) from {1, 2} × {1, 2} × {1, 2}. 1, 1, 1), (1, 1, 2), (2, 2, 1), (2, 2, 2
Iterating over all the constraints of CS will eventually exhibit the inconsistency of this example.

In fact, this shows that exact filtering of a CS over D(x 1 )× ..× D(x n ) can be reached if one computes over a complete lattice built over the set of possible valuations: (P(D(x 1 ) × .. × D(x n )), ⊆, , ). This lattice will be called the concrete lattice in the rest of the paper. Of course, computing over the concrete lattice is usually unreasonable, as it requires to examine every tuple of the Cartesian product D(x 1 ) × .. × D(x n ) w.r.t. consistency of each constraint.

Domain-consistency filtering

For binary constraint systems, the most successful local consistency filtering is arc-consistency, which ensures that every value in the domain of one variable has a support in the domain of the other variable. The standard extension of arc-consistency for constraints of more than two variables is domainconsistency (also called hyper-arc consistency [START_REF] Marriott | Programming with Constraints : An Introduction[END_REF]). Roughly speaking, the abstraction that underpins domain-consistency filtering aims at considering each variable domain separately, instead of considering the Cartesian product of each individual domain. More formally,

Constraint-based reachability

The concretization function is a function γ arc : D ♯ arc → P(D) such that

γ arc ((S 1 , . . . , S n )) = S 1 × . . . × S n
If ⊑ ♯ arc , ♯ arc and ♯ arc denote respectivelly the inclusion, union and intersection of two tuples of sets, then we got the following Galois connection:

(P(D), ⊆, , ) γ arc ⇆ α arc (D ♯ arc , ⊑ ♯ arc , ♯ arc , ♯ arc )
The proof follows comes the monotonicity of the projection and Cartesian product. From Theorem 3, we get:

Definition 7
The best sound approximation of the exact filtering operator f i is

f ♯ i arc def = α arc • f i • γ arc Theorem 7
Let p be a filtering operator associated with constraint c i , then p computes domain-consistency iff p = f ♯ i arc .

This theorem implies that domain-consistency is the strongest property that can be guaranteed by a filtering operator using the abstraction α arc . A proof is given in the Appendix of the paper. Let us consider now the function

f ♯ arc such that f ♯ arc = f ♯ 1 arc • . . . • f ♯ n arc . As f ♯ arc is a sound approximation of f C then sol(C) = gfp( f C ) ⊆ γ arc (gfp( f ♯ arc ))
This result shows if necessary that constraint propagation over domain-consistency filtering operators computes an over-approximation of the solution set of C.

Bound-consistency filtering

Following the same scheme, AI can be used to show the abstraction that underpins constraint propagation with bound-consistency filtering (also called interval-consistency). But, firstly, let us recall the definition of bound-consistency we consider in this paper, as several definitions exist in the literature [START_REF] Chiu | Finite Domain Bounds Consistency Revisited[END_REF] :

Definition 8 (Bound-consistency) A domain D is bound-consistent for a constraint c where vars(c) = {x 1 , .., x n } iff for each variable x i , 1 ≤ i ≤ n and for each d i ∈ {in f D x i , sup D x i } there exist integers d j with in f D x j ≤ d j ≤ sup D x j , 1 ≤ j ≤ n, j = i such that σ = {x 1 → d 1 , .., x n → d n } is an integer solution of c.
Roughly speaking, this approximation considers only the bounds of the domain of each variable and approximates each domain with an interval. Let I (S) = [min(S), max(S)] be the smallest interval that contains all the elements of a finite set of integers S. Similarly, I -1 (I) denotes the set of integers of an interval I :

I -1 ([a, b]) = {x ∈ Z | a ≤ x ≤ b}.
The abstract domain we consider for bound-consistency is D ♯ bound = I (P(D(x 1 )))×. . .×I (P(D(x n ))). Given a tuple of sets (S 1 , . . . , S n ) and a tuple of intervals (I 1 , . . . , I n ), we consider the functions α inter and γ inter such that: α inter (S 1 , . . . , S n ) = (I (S 1 ), . . . , I (S n )) γ inter (I 1 , . . . I n ) = (I -1 (I 1 ), . . . , I -1 (I n ))

Let α bound : P(D) → D ♯ bound be an abstraction function such that α bound = α inter • α arc and γ bound : D ♯ bound → P(D) be a concretization function such that γ bound = γ arc • γ inter If ⊑ ♯ bound , ♯ bound and ♯ bound respectively denote inclusion, union and intersection of intervals (component by component) then we get the following Galois connection:

(P(D), ⊆, , ) γ bound ⇆ α bound (D ♯ bound , ⊑ ♯ bound , ♯ bound , ♯ bound )
Let f ♯ i bound be the most accurate sound approximation of f i , then we get:

f ♯ i bound = α bound • f i • γ bound = α inter • f ♯ i arc • γ inter

Theorem 8 If p is a filtering operator associated to constraint c i , then p computes bound-consistency

iff p = f ♯ i bound .
This theorem, proved in Appendix, implies that bound-consistency is the strongest property that can be reached with an operator based on the α bound abstraction. Consider now the function

f ♯ bound such that f ♯ bound = f ♯ 1 bound • . . . • f ♯ n bound . As f ♯ bound is a sound approxi- mation of f C , then sol(C) = gfp( f C ) ⊆ γ bound (gfp( f ♯ bound ))
This result shows if necessary that constraint propagation based on bound-consistency computes a sound over-approximation of the solution set of C. In addition, as f ♯ bound is also a sound over-approximation of f ♯ arc , then γ arc (gfp( f ♯ arc )) ⊆ γ bound (gfp( f ♯ bound )) meaning that filtering with bound-consistency provides an over-approximation of the results given by a filtering with domain-consistency.

New filtering consistencies based on abstract domains

In the previous section, classical filtering consistencies are interpreted in terms of abstract domain computations. In this section, we propose a new filtering consistency based on the Polyhedral abstract domain [START_REF] Cousot | Automatic Discovery of Linear Restraints Among Variables of a Program[END_REF].

Linear relaxations

When non-linear constraints are involved in a constraint store, approximating them with linear constraints is natural in order to benefit from powerful Linear Programming techniques. These techniques can be used to check the satisfiability of the constraint store when the approximation is sound. If the approximate constraint system is unsatisfiable so is the non-linear constraint system. But, in the context of optimization problems, the approximation can also be used to prune current bounds of the function to optimize.

Another form of approximation comes from the domain in which the computation occurs. A linear problem over integers can be relaxed in the domain of rationals or reals and solved within this domain. As the set of integers belongs to the rationals and reals, an integer solution of the relaxed problem is also a solution of the original integer problem, but the converse is false. In this paper, we will consider both kinds of approximations under the generic term of "linear relaxations".

Computing a linear relaxation of a constraint system CS aims at finding a set of linear constraints that characterizes an over-approximation of the solution set of CS. It is not unique but for trivial reasons, we are more interested in the tighter possible relaxations. The tightest linear relaxation is the convex hull of the solution set of CS but computing this relaxation is as hard as solving CS. For CS over finite domains, the problem is therefore NP hard. Whenever a relaxation is computed by using the current bounds of variable domains, it is called dynamic and the consistencies presented in the rest of the section are compatible with dynamic linear relaxations.

Polyhedral-consistency filtering

Let Poly be the abstract domain of closed convex polyhedra with rational coefficients. As said previously, Poly is not a complete lattice, and then we cannot define a Galois connection between Poly and the lattice of the solutions. Nevertheless, the concretization function γ poly : Poly → P(D) can be defined as the function that returns the integer points of a given polyhedron:

γ poly (S ♯ ) = int sol(S ♯ )
Here, int sol stands for the whole set of integer solutions of a set of linear constraints As S ♯ is bounded, γ poly (S ♯ ) is finite.

Without a Galois connection, we do not expect the polyhedral-consistency proposed in this section to be optimal w.r.t. the abstract domain. Hence, we only show that the filtering algorithm that computes this consistency is a sound approximation of the exact filtering operator. Definition 9 Let α box be the following abstraction function

α box : D ♯ bound → Poly such that α box (([a 1 , b 1 ], . . . , [a m , b m ])) = {a 1 ≤ x 1 ≤ b 1 , . . . , a m ≤ x m ≤ bm}
and the concretization function γ box : Poly → D ♯ bound :

γ box (P) =    ([⌈min(x 1 , P)⌉, ⌊max(x 1 , P)⌋], . . . , [⌈min(x m , P)⌉, ⌊max(x m , P)⌋]) if ∀i, ⌈min(x i , P)⌉ ≤ ⌊max(x i , P)⌋ / 0 otherwise
where ⌊x⌋ (resp. ⌈x⌉) stands for the next smallest (resp. largest) integer of x, and min(v, P) ( resp. max(v, P)) computes the smallest (resp. largest) value of v corresponding to a point of P. Both α box and γ box link the polyhedral abstract domain with the interval abstract domain. The abstraction function α box maps a set of intervals into a polyhedron by adding two inequalities per variable, while the concretization function γ box maps a polyhedron into a set of intervals by computing first the smallest hypercuboid containing the polyhedron and second the greatest hypercuboid with integer bounds. The behaviour of these two functions is illustrated in Fig. 2.

α box γ box
(c) = {x 1 , .., x n } iff for each variable x i , 1 ≤ i ≤ n and for each d i ∈ {in f D x i , sup D x i } there exist rationals r j with in f D x j ≤ r j ≤ sup D x j , 1 ≤ j ≤ n, j = i such that σ = {x 1 → r 1 , .., x n → r n } is a (rational) solution of a linear relaxation of c.
The rationale behind this definition is to benefit from efficient polyhedral techniques over the rationals to filter the variation domain of variables. Of course, interesting implementations of this filtering consistency should trade between efficiency and precision as integer linear constraint solving is costly (NP hard problem) even for bounded domains. It is worth noticing that the definition depends on the quality of the underlying linear relaxation. On the one hand, a linear relaxation which over-approximate c by True (the whole search space) is useless while on the other hand a linear relaxation which exploits piecewise over-approximations of c is often too costly. We give examples of polyhedral-consistency filtering in function of various linear relaxations.

Note that D is bound-consistent for all the constraints of CS.

The simplest linear relaxation that can be considered is the one that ignores non-linear constraints. In this example, c is over-approximated by True and then D viewed as x ≥ -7, x ≤ 10, y ≥ -7, x ≤ 10, z ≥ 3, x ≤ 10, z = x + y is then polyhedral-consistant w.r.t. this linear relaxation. Note that this approach can be generalized by associating a new fresh variable to the non-linear term x * y with a domain computed using the bounds x and y. In this example, this does not help but it could help on other examples. Another linear relaxation consists in building a polyhedron from the "bounds" of x * y in D = x ∈ -7..10, y ∈ -7..10, z ∈ 3..10. By considering the 2-dimensional polyhedron {(1, 10), (10, 1), (-1, -7), (-7, -1)} we get that a linear relaxation of c in domain

D is 11x -8y + 69 ≥ 0 -x -y + 11 ≥ 0 -8x + 11y + 69 ≥ 0 x + y + 8 ≥ 0
Filtering with the polyhedral-consistency, we get that x ∈ -2..9, y ∈ -2..9, z ∈ 3..10 where D(x) and D(y) have been pruned. These results can be easily computed using a Linear Programming tool and truncation operators. For example, using the clpq library of SICStus Prolog which implements a simplex over the rationals, the following request permits to compute the max bound of variable x: {X >= -7, X =< 10, Y >= -7, Y =< 10, Z >= 3, Z =< 10, Z = X+Y, 11 * X -8 * Y+ 69 >=0, -X -Y + 11 >= 0, -8 * X + 11 * Y +69 >= 0, X + Y + 8 >=0}, sup(X, R).

R = 179/19

% then max bound of x is 9

Finally, we can automate the computation of linear relaxations of c by considering the following trivial constraints, which are always true for any x and y:

(x -in f D x)(y -in f D y) ≥ 0 (x -sup D x)(y -in f D y) ≤ 0 (x -in f D x)(y -sup D y) ≤ 0 (x -sup D x)(y -sup D y) ≥ 0
By decomposing these constraints, using the original bounds of x, y, z and replacing the quadratic term x * y by z, we get:

7x + 7y + z + 49 ≥ 0 10x -7y -z + 70 ≥ 0 -7x + 10y -z + 70 ≥ 0 -10x -10y + z + 100 ≥ 0
Filtering with the polyhedral-consistency, we get that x ∈ -2..9, y ∈ -2..9, z ∈ 3..10 where D(x) and D(y) have been pruned. These domains are still bound-consistent but another tighter relaxation can be computed with these new bounds:

2X + 2Y + Z + 4 = 0 9X -2Y -Z + 18 = 0 -2X + 9Y -Z + 18 = 0 -9X -9Y + Z + 81 = 0
and then filtering again permits to get that x ∈ 0..8, y ∈ 0..8, z ∈ 3..10. Here, filtering by bound-consistency leads to prune the domains to: x ∈ 1..8, y ∈ 1..8, z ∈ 3..10. Then, by iterating these two process, we get the only solution to CS which is: x ∈ 2..2, y ∈ 2..2, z ∈ 4..4. This showed how dynamic linear relaxations can be used to solve a non-linear CS.

The w constraint operator

In this section, we present the w constraint operator which captures iterative computations, and how it is processed by a constraint solver. The constraint operator has been introduced a long time ago in [START_REF] Gotlieb | Automatic Test Data Generation Using Constraint Solving Techniques[END_REF][START_REF] Gotlieb | A CLP Framework for Computing Structural Test Data[END_REF] and was further refined using Abstract Interpretation (AI) techniques [START_REF] Denmat | An Abstract Interpretation Based Combinator for Modeling While Loops in Constraint Programming[END_REF]. In the following, we recall its interface and semantics and show how fixed point computations can be used to filter inconsistant values of the underlying relation. We also explain how the Polyhedral abstract domain is used to approximate the fixed point computations.

w as a relation over memory states

The w operator captures a relation over three memory states that represent the state before, within and after the execution of an iterating statement. In this paper, we do not specify what a memory state is, or what the iterating statement is, as the approach is generic regarding the content of a memory state and the concrete syntax of the iterator. However, in order to ease the understanding, the reader can consider a memory state to be a mapping between variables of the program to values. More complex examples of memory states in relation with w can be found in [START_REF] Charreteur | Modelling dynamic memory management in Constraint-Based Testing[END_REF] and [START_REF] Charreteur | Constraint-Based Test Input Generation for Java Bytecode[END_REF]. The relation w is expressed with the following syntax: w(M 1 , M 2 , M 3 , Dec, Body)} where M 1 denotes the memory state before execution of the iteration, M 2 denotes the memory state reached at the end of execution of the Body, while M 3 denotes the state after execution, Dec is a boolean syntactical expression, and Body is a list of statements. This three-states consideration is inspired by the Static Single Assignment of a program [START_REF] Wegman | Constant Propagation with Conditional Branches[END_REF]. If the state of M is irrelevant for a given computation, we simply write . Note that Body may also contain other iterators, and thus w is meant to be a compositional operator.

The semantics of w is the semantics of an iterating statement (i.e., repetitive application of Body over an input state, while Dec is true).

We note w n = n w • w . . . w where • is the application composition.

Background on w

As described in [START_REF] Gotlieb | A CLP Framework for Computing Structural Test Data[END_REF], the operational semantics of w within a constraint solver is expressed as a set of guarded-constraints:

{(C 1 -→ C 2 ) i } 1≤i≤n .
If C 1 is entailed by the constraint store then C 2 is added to it, and the relation w is solved. If C 1 is disentailed, then the guarded-constraint is discarded and no more considered in further analysis. Finally, if none of these (dis-) entailment deductions is possible, the guarded-constraint just suspends in the constraint store. The set of guarded-constraints is considered each time the constraint w awakes in the constraint store, so that it captures the essence of the iteration through rewriting in recursive calls. In addition, substitution of variables must be considered to faithfully represent the constraints in a w relation. Dec M 3 ←M 1 simply denotes the constraint Dec where program variables from M 3 have been substituted by the variables from M 1 . With these notations, the w relation is expressed as follows: w(Dec, M 1 , M 2 , M 3 , Body) iff

• Dec M 3 ←M 1 -→ Body M 3 ←M 1 ∧ w(Dec, M 2 , M new , M 3 , Body M 2 ←M new ) • ¬(Dec M 3 ←M 1 ) -→ M 3 = M 1 • ¬(Dec M 3 ←M 1 ∧ Body M 3 ←M 1 ) -→ ¬(Dec M 3 ←M 1 ) ∧ M 3 = M 1 • ¬(¬Dec M 3 ←M 1 ∧ M 3 = M 1 ) -→ Dec M 3 ←M 1 ∧ Body M 3 ←M 1 ∧ w(Dec, M 2 , M new , M 3 , Body M 2 ←M new ) • join(Dec M 3 ←M 1 ∧Body M 3 ←M 1 ∧w(Dec, M 2 , M new , M 3 , Body M 2 ←M new ), ¬(Dec M 3 ←M 1 )∧M 3 = M 1 )
The two former guarded-constraints implement forward analysis, by examining the entailment of Dec. Depending on the entailment of Dec, a recursive call to a new w is added to the constraint store. The two followings implement backward reasoning by examining the differences between the stores after and before execution of the iteration. Finally, the last operation, called join, is the most tricky one and implements union of stores in case of suspension of the operator. This join operation is realized iff none of the previous guarded-constraints has been solved. The rest of the Section is devoted to the presentation of this operator, which is implemented as an abstract operation over abstract domains.

Concrete fixed point computation

For a given w operator, let T be the following set:

T = {(M i , M j ) | ∃k | w k (M i , M j , , Dec, Body)}
T represents all pairs of memory states that are in relation through the w statement, but still, not all those pairs can be considered as solutions of the relation, as some pairs can only be reached in temporary states of the execution. For this reason, we introduce the set Z w :

Z w = {(M i , M j ) | (M i , M j ) ∈ T ∧ M j ∈ sol(¬Dec)}
where sol(C) denotes the set of solutions of a constraint C. T can be seen as the least fixed point of:

T i+1 = {(M k , M j ) | (T i ∧ w(M k , , M j , Dec, Body))} ∪ T i (4) T 0 = {(M 1 , M 1 )} (5)
and Z w can be computed by filtering the pairs of the fixed point.

For instance, considering

M 1 = x → 0 ∨ x → 1 ∨ x → 2 ∨ x → 3 and w(M 1 , M 2 , M 3 , x < 2, x = x + 1)
, and using the notation (0, 0) for denotating (x → 0, x → 0), the fix point computation is as follows:

T 0 = {(0, 0), (1, 1), (2, 2), (3, 3)} T 1 = {(0, 1), (1, 2)} ∪ T 0 = {(0, 0), (0, 1), (1, 1), (1, 2), (2, 2), (3, 3)} T 2 = {(0, 1), (0, 2), (1, 2)} ∪ T 1 = {(0, 0), (0, 1), (0, 2), (1, 1), (1, 2), (2, 2), (3, 3)}

T 3 = T 2
Consequently, the solutions set Z w of w(M 1 , M 2 , M 3 , x < 2, x = x + 1) is:

Z w = {(a, b) | (a, b) ∈ T 3 ∧ (x → b) ∈ sol(x ≥ 2)} = {(0, 2), (1, 2), (2, 2), (3, 3)}
Computing Z w is undecidable in general as there is no termination guarantee of the iterating process. This is the reason why this computation is usually abstracted using abstract domain computation.

Abstracting the fixed point computation

Implementing the join operator mentionned above can be done by abstracting the computation of the fixed point within the Polyhedral abstract domain. Let P ♯ be a conjunction of linear restraints, the intersection of which defines a convex polyhedron, that over-approximates the set T . Hence, we can compute P ♯ as the least fixed point of:

P i+1 = {(M k , M j ) | (P i ∧ α poly ( w(M k , M j , , Dec, Body))) ⊔ P i ( 6
)
P 0 = {(α poly ((M 1 , M 1 )) (7)
Compared to eq. 4 and 5, the computation is realized in the abstract domain using α poly the abstraction function of the Polyhedral abstract domain.

Let Z ♯ w be the approximation of the set of solutions of w, obtained by application of α poly :

Z ♯ w = {(M i , M j ) | (M i , M j ) ∈ P ♯ ∧ M j ∈ α poly (sol(¬Dec))}
Looking at the above example where M is just composed of the mapping of x → v, it is worth introducing different representations of the stores as we progress in the fixed point computation. When P i is computed M k and M j that contains x j , we note: P i (x k , x j ). If P i is then considered over y k , y j , then we will simply write P i (y k , y j ) and apply variable substitution. With these notations, we have the following computation:

P 0 (x in , x out ) = x in ≥ 0 ∧ x in ≤ 3 ∧ x in = x out P 1 (x in , x out ) = (P 0 (x in , x 0 ) ∧ x 0 ≤ 1 ∧ x out = x 0 + 1) x in ,x out ⊔ P 0 (x in , x out ) = (x in ≥ 0 ∧ x in ≤ 1 ∧ x out = x in + 1) ⊔ P 0 (x in , x out ) = x in ≥ 0 ∧ x in ≤ 3 ∧ x out ≤ x in + 1 ∧ x out ≥ x in P 2 (x in , x out ) = (P 1 (x in , x 1 ) ∧ x 1 ≤ 1 ∧ x out = x 1 + 1) x in ,x out ⊔ P 1 (x in , x out ) = (x in ≥ 0 ∧ x in ≤ 3 ∧ x in ≤ x out -1) ⊔ P 1 (x in , x out ) = x in ≥ 0 ∧ x in ≤ 3 ∧ x out ≤ x in + 2 ∧ x out ≥ x in ∧ x out ≤ 4 P 3 (x in , x out ) = (P 2 (x in , x 2 ) ∧ x 2 ≤ 1 ∧ x out = x 2 + 1) x in ,x out ⊔ P 2 (x in , x out ) = (x in ≥ 0 ∧ x in ≤ 3 ∧ x in ≤ x out -1) ⊔ P 2 (x in , x out ) = P 2 (x in , x out )
Fig. 3 illustrates the difference between the abstract fixed point and the approximate fixed point. Points in the figure correspond to the elements of T 3 , while the grey zone represents the convex polyhedron defined by P 3 . An approximation of the solutions of w(M 1 , M 2 , M 3 , x < 2, x = x + 1) is given by:

Q = P 3 (x 1 , x 3 ) ∧ x 3 ≥ 2 = x 3 ≥ 2 ∧ x 3 ≤ 4 ∧ x 1 ≤ x 3 ∧ x 1 ≤ 3 ∧ x 1 ≥ x 3 - 2
On the Polyhedral domain, convergence of the fixed point computation over w(M 1 , M 2 , M 3 , x < 2, x = x + 1) can be enforced by using widening techniques. The computation of P k+1 is modified in order to use a widening operator ∇ [START_REF] Cousot | Automatic Discovery of Linear Restraints Among Variables of a Program[END_REF]. Thus, we have:

P k+1 = P k (Init, Out)∇(P k ∧ α poly (w(M 1 , M 2 , M 3 , Dec, Body)))
A concrete algorithm for computing this approximation is given in [START_REF] Denmat | An Abstract Interpretation Based Combinator for Modeling While Loops in Constraint Programming[END_REF], which permits to build implementation of w in a constraint solver. As rooted in the Abstract Interpretation domain, the relation w inherits from some of its fundamental correctness results, i.e., soundness and termination. However, it is worth pinpointing some differences.

Usually, a convex abstract polyhedron denotes the set of linear relations that hold over variables at a given point of a sequential program under analysis. As the goal here is to correctly approximate the set of solutions of a w relation, the polyhedron describes relations between input and output values and, thus, they involve more variables in the equations. In Abstract Interpretation, the analysis can be performed only once, whereas, in the case of the w relation, the join operation is launched everytime the relation is awaked without being succesfull in solving one of the guarded-constraint. As a consequence, we found out that it was not reasonable to use standard libraries to compute over polyhedra, such as PPL [START_REF] Bagnara | Possibly Not Closed Convex Polyhedra and the Parma Polyhedra Library[END_REF], because they use a dual representation for Polyhedra, which is a source of exponential time computations for the conversion.

Illustrative example

Looking at an iterative computation over unbounded domains as a relation captured by a w constraint operator is interesting for adressing Constraint-Based Reacheability problems. On the one hand, the suspension mechanism offered by constraint reasoning allows us to cope with the approximation problem, i.e., the set of states that is considered is determined by the informations existing in the constraint store, which makes the reasoning more accurate w.r.t. the property to be demonstrated. On the other hand, adding abstract domain computations to the w relation allows us to increase the level of deductions that can be achieved at each awakening of the w constraint operator. To illustrate this remark, consider the following C program: f( int i, ... ) { a. j = 100; b.

while( i > 0) c.

{ j=j+1 ; i=i-1 ;} d.

... e. if( j > 500) f. ...

A typical reachability problem is to find out a value of i such that statement f. is executed. Existing approaches for solving this reachability problem consider a path passing through f., e.g., a-b-d-e-f, and try to solve the path condition attached to this path. In this case, it means extracting constraint j 1 = 100 ∧ i 1 ≤ 0 ∧ j 1 > 500 and solving it to show that the constraint system is unsatisfiable, i.e., the corresponding path is infeasible. Then, these approaches backtrack to select another path (e.g., a-b-c-b-d-e-f with path condition

j 1 = 100 ∧ i 1 > 0 ∧ j 2 = j 1 + 1 ∧ i 2 = i 1 -1 ∧ i 2 ≤ 0 ∧ j 2 > 500)
and repeat the process again, until a satisfiable path condition is found. This example is pathologic for these approaches, as only the paths that iterate more than 400 times in the loop will reach statement f.. Hopefully, using the constraint operator w(M 1 , M 2 , M 3 , i > 0, j = j + 1 ∧ i = i -1) permits us to unrool dynamically 400 times the loop without backtracking. The relational analysis performed on the Polyhedral abstract domain by the w operator determines that j outi in = 100 whatever be the number of loop unrollings. Here, combining precise constraint reasoning in the concrete domain, with constraint extrapolation through abstract domain computations, offers us an efficient way of solving reachability problems on infinite-state systems.

Conclusions

In this paper, we have presented Constraint-Based Reachability as a process to combine constraint reasoning and abstraction techniques for solving reachability problems in infinite-state systems. The contribution is two-fold: first, we have revisited constraint consistency-filtering techniques by the prism of abstract domain computations ; second, we explained how to introduce abstract domain computation within the w constraint operator reasoning. We have illustrated these notions with several examples in order to ease the understanding of the reader. This appraoch has been implemented and tested on several problems, including real-world programs [START_REF] Gotlieb | EUCLIDE: A Constraint-Based Testing platform for critical C programs[END_REF][START_REF] Gotlieb | TCAS software verification using Constraint Programming[END_REF]. The goal is now to broader the scope of these techniques that combine constraint reasoning and abstraction techniques, to adress fundamental problems such as reachability in infinite-state systems.

Definition 3 (Figure 1 :

 31 Figure 1: Static approximations of fixed point computations in complete lattices

Theorem 3 (

 3 Smallest sound approximation) Let (D, ⊑, , ) γ ⇆ α (D ♯ , ⊑ ♯ , ♯ , ♯ ) be a Galois connection, and a function f : D → D, then the smallest sound approximation of f is α • f • γ

Definition 5 (

 5 Sound approximation without a Galois connection) Let (D, ⊑) and (D ♯ , ⊑ ♯ ) be two posets, let γ : D ♯ → D be a monotonic function and f : D → D a function, then the function f

Definition 6 (

 6 Domain-consistency) A domain D is domain-consistent for a constraint c where vars(c) = {x 1 , .., x n } iff for each variable x i , 1 ≤ i ≤ n and for each d i ∈ D(x i ) there exist integers d j with d j ∈ D(x j ), 1 ≤ j ≤ n, j = i such that σ = {x 1 → d 1 , .., x n → d n } is an integer solution of c.Consider the domains D = D(x 1 ) × .. × D(x n ) and D ♯ arc = P(D(x 1 )) × . . . × P(D(x n )) and the abstraction function α arc : P(D) → D ♯ arc which maps S ∈ P(D) to α arc (S) = ({x 1 | x ∈ S}, . . . , {x n | x ∈ S})

Figure 2 :

 2 Figure 2: Connection between the Polyhedral and Interval abstract domains

Example 2

 2 Consider the following CS: z = x + y, z = x * y, let c be the second constraint of CS: c = (z = x * y) and let D be x ∈ -7..10, y ∈ -7..10, z ∈ 3..10.

Figure 3 :

 3 Figure 3: Exact and approximated fixed point

Program semantics captures formally all the possible behaviours of a program.

The union of two polyhedra is not a polyhedron, this is the reason why convex hull or any relaxation of it must be employed.

Acknowledgements

We are indebted to Bernard Botella and Mireille Ducassé for fruitful discussions on earlier versions of this work.

Appendix

This appendix contains the proofs of some of the results stated in the paper.

Theorem 9 Let p be a filtering operator associated with constraint c i , then p computes domain-consistency

. From the definitions of f i and γ, we get that S 1 is the solution set of constraint c i , given the initial domains S (we write S 1 = sol(c i , S)). Hence, S ′ = α arc (S 1 ) = (A 1 , . . . , A m ) with A k = {x k | x ∈ sol(c i , S)}. So, c i computes domain-consistency. (⇒) Let p be a domain-consistency filtering operator. Suppose that there exists S such that p(S) = (A 1 , . . . , A m ) be strictly greater than f ♯ i arc (S) = (B 1 , . . . , B m ). Then, there exists at least one k such as A k B k . Hence, there exists an element x k of A k that does not belong to any solution of constraint c i . Hence, p cannot computes domain-consistency which is contradictory with the hypothesis. On the other side, p cannot be smaller than f ♯ i arc as it means that the filtering operator removes solutions. Hence, if p computes domain-consistency then p = f ♯ i arc . 2 Theorem 10 If p is a filtering operator associated to constraint c i , then p computes bound-consistency iff p = f ♯ i bound . Proof 2 (⇐) From theorem 9, given initial intervals I, the domains f ♯ i arc •γ inter (I) are domain-consistent for constraint c i . Applying function α inter is similar to the process that keeps extremal values of each ele- ment of f ♯ i arc • γ inter (I). Hence, the resulting intervals satisfy the bound-consistency property. (⇒) (similar to the proof of theorem 9) If the filtering operator p is greater than f ♯ i bound , then the computed intervals contain at least one bound that is not part of a solution of c i , violating so the boundconsistency property. On the contrary, by supposing that p is smaller than f ♯ i bound then solutions are lost and p is no more a filtering operator. Hence, if p is a filtering operator guaranteeing bound-consistency then p = f ♯ i bound . 2