Vincent Dahirel

Xudong Zhao

Baptiste Couet

Guillaume Batôt

Marie Jardat

Hydrodynamic interactions between solutes in multiparticle collision dynamics

Multiparticle collision dynamics (MPCD) enables to simulate fluid dynamics including both hydrodynamics and thermal fluctuations. Its main use concerns complex fluids, where the solvent interacts with concentrated solutes, may they be colloidal particles, polymers or electrolytes. A key difficulty concerns the way one couples the fluid to the solute particles, without losing the key advantages of the MPCD method in therm of computational efficiency. In this paper, we investigate the dynamical properties of solutes that are coupled to the fluid within the collision step, i.e. when local momentum exchange between fluid particles occur. We quantify how the volume where momentum exchange is performed (the size of the collision cells) constrains the hydrodynamic size of the solute. Moreover, we show that this volume should be taken smaller than the structural size of the solutes. Within these constraints, we find that the hydrodynamic properties of an electrolyte solution are similar to the behaviour predicted by the Fuoss-Onsager theory of electrolyte dynamics. However, it is also clear that mapping the diffusion time scale to that of a real system cannot be done quantitatively with this methodology.

I. INTRODUCTION

The dynamics of solutes in a solution drives many properties, from their thermic or electrical conductivity to the rate of chemical reactions. Dynamical properties of solutes depend on hydrodynamic interactions mediated by the solvent. These interactions play a role in a very large concentration range. Typically, in electrolytes 1 , they become significant at concentrations larger than 10 -2 mol.L -1 , and in some case they become screened in very dense solutions, such as polymer melts 2 . Most transport theories of simple solutions, electrolyte solutions, dilute polymers, or colloidal suspensions rely on a description of hydrodynamic interactions derived from the Stokes equation of fluid dynamics at low Reynolds numbers. The Fuoss-Onsager theory of electrolyte transport 1,3 and the classical theories of polymer transport 2,4-6 both include hydrodynamics using the Oseen tensor. In this framework, the velocity of a solute is influenced by the velocities of surrounding solutes through a tensor depending on the distances between particles. The Oseen tensor can be derived by evaluating the effect of a force applied on a fixed point of the fluid described by the Stokes equation. The size and shape of the particles have no explicit influence on hydrodynamic interactions within this modelling. Finding numerical alternatives to these limited theoretical treatments without resorting to atomistic numerical simulations has given rise to a variety of mesoscopic methods. Brownian Dynamics with the Rotne-Prager hydrodynamic tensor 6 is the simulation technique that is the most natural 8 , as it contains ingredients from classical theoretical treatments. This method has been very successful in computing transport coefficients for various systems [9][10][11][12][13][14][15][16] but it has several major pitfalls. First, when the system is too concentrated, the Rotne-Prager tensor can be non positive definite for some configurations of the system, and the simulation cannot proceed. This is a clear limitation if one is interested in crowded media, such as the interior of a biological cell. Secondly, when strong attractive interactions between particles exist, random displacements can lead the system in regions of phase space that should not be realistically explored, which leads to strong instabilities of the simulation. This difficulty can be partially overcome by the use of Metropolis Adjusted Langevin Algorithm (MALA), see e. g. ref. 17 . Moreover, it is very difficult to adapt such simulation strategies to concentrated solutions or confined solutions for which the use of the Rotne-Prager tensor is no longer valid 7 .

Two popular alternatives to Brownian Dynamics are Dissipative Particle Dynamics 18 , and Multi-Particle Collision Dynamics (MPCD) 19,[START_REF] Gompper | In Advanced Computer Simulation Approaches for Soft Matter Sciences III[END_REF] . In the present paper, we focus on the use of the latter to study systems with strong hydrodynamic interactions between solutes. In MPCD, an explicit but highly simplified description of the solvent is used, where ballistic motions and local momentum exchanges between solvent particles are tuned to reproduce the properties of a fluid at the Navier-Stokes level. It is a discrete solvent with the hydrodynamic properties of a continuous solvent, whose hydrodynamic regime can be chosen thanks to a relatively small number of parameters. Solute particles can be embedded in the MPCD solvent bath, and evolve through a classical molecular dynamics algorithm; Hydrodynamic interactions between solutes emerge in this case. One advantage of MPCD is that it describes hydrodynamics more adequately than Brownian Dynamics with the Rotne-Prager tensor, and that it can in principle be used for any boundary conditions. Moreover, the algorithm is particularly suited to parallelization as the local momentum exchanges happen in collision cells defined by a grid of fixed size a 0 .

There are several ways to couple the solute particles to the MPCD fluid. A central repulsive force between the solute and the solvent can be added. It creates a spherical zone around the solute depleted from the solvent [START_REF] Malevanets | [END_REF]22 . This scheme leads to an effective slip boundary condi-tion at the surface of the solute. Alternative schemes can lead to an effective stick boundary condition [23][24][25] . In cases where the solute particles are rather small, another possibility is to couple the solute with the solvent during the collision step, when the momentum exchange occurs. This coupling scheme is rather crude, but it is the most efficient from the computational point of view. In the following, we refer to this coupling scheme as collisional coupling (denoted by MPCD-CC in what follows). This scheme has been used for instance by Winkler and coworkers to study small polymer chains 26,27 , and allowed them to recover the scaling laws of Zimm theories (see also Ref. 28 ). More recently, this method has been used for a model of solutes with several sites 29 , in order to mimic hydrodynamic interactions at the surface of a colloidal particle. Despite these successes, some of us found in a previous work that the diffusion coefficients of simple electrolytes predicted by MPCD-CC differed from those obtained by Brownian Dynamics with hydrodynamic interactions, in contrast with more expensive coupling algorithms 30 . Indeed, one limitation of the MPCD-CC scheme is that the solute influences the solvent at the length scale of the collision cell a 0 only, so that the hydrodynamic size of the solute is of the order of a 0 in any case. In real systems, the hydrodynamic radius of a simple ion or of a nanoparticle is close to its structural radius derived from the atomic structure. It means that the structural size should also be of the order of a 0 , or in other word that the size of the cell in MPCD should be chosen to match the structural size when the collisional coupling is used. In practice, in several papers, authors suggest to choose a value of the structural diameter of the order of a 0 [31][32][33] , but there is no quantification of the effect of the cell size on the transport properties.

In the present article, we focus on the MPCD-CC method and investigate its ability to reproduce hydrodynamic interactions as well as transport theories at the Oseen level. First, we compute from different methods the effective hydrodynamic radius at infinite dilution of a solute in collisional coupling with the MPCD solvent. We show that it is in any case of the order of one third of the cell size a 0 . If the structural radius of the solute is chosen equal to the hydrodynamic radius, it means that several solutes may be in the same collision cell at a given step, which might lead to spurious effects. We thus also investigate the influence of the cell size on the diffusion coefficient of solutes. The results are compared to the ones we previously obtained for the same systems without hydrodynamic interactions and also with MPCD and a central force between solvent and solute particles. For hard spheres, we find an upper limit for the cell size compared to the structural radius over which the diffusion coefficient spuriously decreases. Finally, we consider the case of solutes with an attractive interaction, as it may increase the probability to find two solutes in the same collision cell. More precisely, we study the transport properties of simple electrolyte solutions. We compare the diffusion coefficient computed by MPCD-CC for two different resolutions of the grid to our previous numerical results, taken here as references. We also compute the electrical conductivity of the solution because i) the electrical conductivity of electrolytes is known to be strongly affected by hydrodynamic couplings, ii) a reliable semi-analytical transport theory accounting for hydrodynamic couplings at the Oseen level exists, able to predict the electrical conductivity of simple electrolytes over a wide range of concentration. We find an excellent agreement between the electrical conductivity computed from MPCD-CC and the transport theory, which shows that in the range of parameters chosen, MPCD-CC is able to capture the hydrodynamic interactions between ions. However, we also show that the constraints on the choice of parameters in the MPCD-CC method prevent us from representing any given real system.

The paper is organized as follows. In Section II we shortly describe the simulation methods and the semianalytical theory used to compute the electrical conductivity. We compute the effective hydrodynamic radius of a solute in collisional coupling with a MPCD solvent in Sec. III. Then, the influence on diffusion coefficients of the size of the collision cell compared to the structural radius of the solute is investigated in Sec. IV for solutions of neutral solutes and for electrolyte solutions. Finally, we compare the electrical conductivity computed by MPCD-CC to the theoretical prediction in Sec. V. The paper ends with a conclusive discussion.

II. METHODS

A. Multi-particle collision dynamics simulations

The fluid in MPCD is represented by pointlike particles, whose positions and velocities evolve in two steps 22 . In the streaming step, positions and velocities are propagated by integrating Newton's equations of motion. Without external forces, this yields a ballistic motion for each fluid particle i:

r i (t + δt c ) = r i (t) + v i (t)δt c ( 1
)
where r i , v i are respectively the position and the velocity of particle i, and δt c is the time step. A second step, the collision step, enables local momentum exchanges between the fluid particles. The simulation box is partitioned into cubic cells of given size a 0 . A randomly oriented axis is defined for each collision cell, and the velocities of fluid particles relative to the velocity of the center of mass of the cell are rotated by an angle α around this axis:

v i (t + δt c ) = v cell c.o.m (t) + R α [v i (t) -v cell c.o.m (t)] (2)
where R α is the rotation matrix and v cell c.o.m the velocity of the center of mass of the cell. The angle α is a fixed parameter. A random shift of the collision grid is performed at each collision step to ensure galilean invariance 34,35 .

The transport properties of the fluid depend on a few parameters: The number of solvent particles per cell γ, the rotation angle α, the time between two collisions δt c [START_REF] Gompper | In Advanced Computer Simulation Approaches for Soft Matter Sciences III[END_REF]22 . It is convenient to use the fluid particle mass m f as the mass unit, the size of the collision cells a 0 as the length unit, and k B T as the energy unit with T the temperature and k B Boltzmann constant. The time unit is then

t 0 = a 0 m f k B T . ( 3
)
The kinematic viscosity of the pure MPCD fluid is ν = ν coll + ν kin 36,37 with

ν coll = 1 λ (1 -cosα) 18 1 - 1 γ + e -γ γ ( 4
)
ν kin = λ 1 [4 -2cosα -2cos(2α)] 5γ (γ -1 + e -γ ) - 1 2 . ( 5
)
where λ is mean free path of fluid particles, λ = δt c /t 0 .

By default, in what follows, the parameters of the MPCD simulations are: {α = 130 • , γ = 5, δt c = 0.1t 0 }. This ensures that the Schmidt number (ratio of the time scale of diffusive mass transfer over the time scale of momentum transfer in the fluid), corresponds to a liquidlike behavior. For this choice of parameters, the kinematic viscosity of the fluid is ν = 0.809 a 2 0 .t -1 0 (with ν coll >> ν kin ), so that the dynamic viscosity η is equal to 4.045 m f a -1 0 .t -1 0 . Solute particles can be immersed in this solvent bath. Their dynamics is then coupled to that of the fluid particles. In what follows, we focus on the so-called collisional coupling scheme, where solute particles i) interact with each other through a given force field, ii) participate to the collision step with solvent particles located in the same cell. We refer to this method as MPCD-CC hereafter. Solutes usually have a mass greater than that of fluid particles. During the streaming step, the position R j and velocity V j of solute j are propagated with the velocity Verlet algorithm often used in standard Molecular Dynamics (MD) simulations:

R j (t + δt M D ) = R j (t) + V j (t)δt M D + F j (t) 2M δt 2 M D , ( 6
)
V j (t + δt M D ) = V j (t) + F j (t) + F j (t + δt M D ) 2M δt M D , ( 7
)
where M is the mass of the solute particle, F j is the force acting on solute j at the beginning of the step, and δt M D is the time step. During the collision step, the velocities of fluid and solute particles are updated following Eq. ( 2) in each collision cell. More details about this simulation scheme can be found in Refs. [START_REF] Gompper | In Advanced Computer Simulation Approaches for Soft Matter Sciences III[END_REF]38 .

Another possible coupling between the dynamics of solutes and the solvent bath consists in adding an explicit short ranged interaction between solutes and fluid particles, preventing fluid particles to penetrate into solutes. In this case, solutes do not participate to the collision step. Details on this simulation scheme, called MPCD-CFC (CFC for Central Force Coupling) in what follows can be found in Ref. 22 . In the present paper, the results obtained with MPCD-CC are compared in a few cases with those obtained in a previous study with MPCD-CFC. For these MPCD-CFC simulations, a purely repulsive short-ranged interaction potential between solutes and solvent was used 30 .

B. Computation of transport coefficients and of the hydrodynamic radius from MPCD-CC

The self-diffusion coefficient of solutes D s is computed using equilibrium trajectories. The mean-squared displacement as function of time is computed, and the diffusion coefficient of solutes is deduced from the slope at long time:

D s = lim t→∞ 1 6t |R j (t + t 0 ) -R j (t 0 )| 2 t0,j . ( 8
)
At infinite dilution, the hydrodynamic radius a hyd can be estimated from the self-diffusion coefficient using the Stokes law, written here in the case of stick boundary conditions between solute and solvent:

D s = k B T 6πηa hyd , (9)
with η the dynamic viscosity of the fluid, η = γν.

The electrical conductivity σ of a solution of charged solutes is computed from equilibrium MPCD simulations using Kubo's formula: (10) with V the volume of the simulation box, N ed the total number of solutes, q i the charge of solute i and V i its velocity.

σ = 1 3k B T V ∞ 0 dt N ed i=1 q i V i (t 0 ) • N ed j=1 q j V j (t 0 + t) t0 ,
We have also used non-equilibrium MPCD simulations to compute the hydrodynamic radius of a solute in collisional coupling. In this procedure, the solute particle is fixed at the center of the simulation box, and we impose a solvent flow along the x-direction. To induce the flow, a velocity with Gaussian distribution centered on a given value v 0 is added to each solvent particle situated in the layer of the simulation box perpendicular to the x-direction, of thickness a 0 , at x = 0. We have chosen v 0 = 0.02 a 0 .t -1 0 . At stationary state, we obtain a flow around the solute with a cylindrical symmetry. In the approximation of laminar flow (Re 1), the analytic solution of the Stokes equation for the fluid around an isolated sphere of radius a hyd fixed at the origin with stick boundary conditions is given by

v r (r, x) = -v ∞ a hyd xr 2(r 2 + x 2 ) 3 2 (11) v θ = 0 v x (r, x) = v ∞ 1 -a hyd 1 2 √ r 2 + x 2 + x 2 2 (r 2 + x 2 ) 3 2
where cylindrical coordinates (r, θ, x) are used, x being the direction of the flow and r the distance to the origin.

The velocity field hence depends on two parameters: the hydrodynamic radius a hyd of the sphere and the fluid velocity far from the sphere v ∞ . The incompressibility of the fluid imposes v ∞ to be equal to the average velocity within a slice of fluid in the plane perpendicular to the flow direction x, whatever the value of x. The value of v ∞ in the flow simulated by MPCD can thus be determined. Then, we have fitted the value of a hyd by minimizing the mean square deviation between the Stokes velocity field (Eq.( 11)) and the velocity field computed by MPCD. As the simulation box has a finite size, the simulated flow field could actually be affected by periodic boundary conditions. Nevertheless, this effect is to the first order proportional to a hyd /L box with L box the size of the simulation box. As we proceed to show, under the conditions of our simulations, a hyd /L box is less than 0.01, so that the finite size effect can be safely neglected.

C. Theoretical computation of the electrical conductivity

The electrical conductivity of electrolyte solutions can be predicted from analytical or semi-analytical theories related to the Fuoss-Onsager theory developed at the beginning of the twentieth century 3 . These theories are based on the assumption that two main effects influence the dynamics of ions in solution, namely electrostatic interactions and hydrodynamic couplings. The forces induce modifications of the velocity of solutes which depend on the structural organisation of the solution through the direct correlation function between two solutes. Analytical expressions of the electrical conductivity have been proposed for the primitive model of electrolytes using the mean spherical approximation as closure relation to solve the Ornstein-Zernike integral equation [START_REF] Hansen | theory of Simple Liquids[END_REF] and to compute the structural correlation functions 1 . In the framework of the primitive model, ions are charged hard spheres embedded in a continuous solvent of given viscosity and dielectric constant. The calculated electrical conductivity was found in excellent agreement with experimental data up to high concentrations 1,[START_REF] Durand-Vidal | [END_REF][41][42] for several electrolytes. We propose here to use this theoretical framework to predict the electrical conductivity of solutions of charged species of given structural and hydrodynamic radii. These theoretical results will be used as references to compare our simulation results with.

More precisely, the velocity of ion i at stationary state under the presence of an electric field E reads

v i = D 0 s k B T q i E + δv hyd i + δv elec i , ( 12
)
where D 0 s is the self-diffusion coefficient of species i at infinite dilution, related to its hydrodynamic radius through the Stokes relation (eq. ( 9)), δv hyd i is the hydrodynamic velocity correction due to hydrodynamic interactions, and δv elec i is the velocity correction due to electrostatic couplings. The hydrodynamic velocity correction δv hyd i can be approximated using the Oseen tensor as 41 :

δv hyd i = j n j q j E 2 3η ∞ 0 rh ij (r)dr ( 13
)
where n j is the density of species i and h ij (r) is the total pair correlation function between species i and j at a distance r. The velocity correction due to electrostatic couplings, δv elec i also depends on the equilibrium pair correlation function, and scales with 1/η. The full expression of this term can be found e. g. in 1,43 . Here, we use the HyperNetted Chain closure relation [START_REF] Hansen | theory of Simple Liquids[END_REF] to compute the total pair correlation function between ions h ij (r). We have checked that for the systems investigated here these correlation functions were very close to those obtained by MPCD. The advantage of the HNC solution for the correlation function is the absence of statistical noise. Moreover, the hydrodynamic radius a hyd and the structural radius a HS are two independent parameters of the calculation, a hyd being related to D • s and a HS being a parameter of the interaction potential. This will allow us to use exactly the same parameters as those of the MPCD-CC simulations to compute the electrical conductivity.

III. EFFECTIVE HYDRODYNAMIC RADIUS OF A SOLUTE IN THE MPCD-CC SCHEME

The effective hydrodynamic radius of a solute can be deduced from its self-diffusion coefficient at infinite dilution from the Stokes law. To compute the diffusion coefficient of a solute in collisional coupling with the MPCD solvent at infinite dilution, we performed simulations at low density of solutes (N ed /L 3 box between 0.0008 for the largest simulation box, and 0.01 for the smallest one, with N ed the number of solutes), without any direct interactions between solutes (the solute-solute correlations are those of a perfect gas). Such simulations require very long trajectories to get enough statistics. The simulations were run for embedded particles of mass M = 10m f , for several box lengths between 10 a 0 and 50 a 0 . Seven independent trajectories of 2.5 × 10 7 steps each, with δt M D = 0.01 t 0 were run. The results are shown in Fig. 1 We observe in Fig. 1 that the self-diffusion coefficient is a linear function of the inverse of the box length due to finite size effects 44,45 . This scaling was first predicted by Hasimoto 44 who studied the velocity flow of a fluid across a periodic cubic array of spheres with L box being the distance between two spheres. For a fluid of viscosity η we expect

D s (L box ) = D s (∞) -2.837 k B T 6πηL box (14) = D s (∞) - 0.0372 L box , ( 15
)
with η = 4.045 m f a -1 0 t -1 0 and k B T = 1. This prediction is in excellent agreement with the fit of our data. It shows that our simulations do correspond to infinite dilution conditions, with a very weak effect of solutefluid interactions on the fluid viscosity, and hence no visible hydrodynamic interactions between solutes inside the simulation box. This also enables us to extrapolate our results to infinite box size so that we obtain D s (∞) = D • s = 0.0422 a 2 0 t -1 0 , and the corresponding hydrodynamic radius of the solute from the Stokes law (eq.( 9)): a hyd = 0.311 a 0 . From this value of the radius, one can compute an effective packing fraction. It is always lower than 4.10 -4 , which confirms that we are in the high dilution regime for indirect hydrodynamic interactions between solutes in the simulation box.

The hydrodynamic radius of a solute in collisional coupling deduced from the diffusion coefficient at infinite dilution is thus found smaller than half the size of the collision cell. One may wonder if we could vary the effective FIG. 2: Left: Self-diffusion coefficient at infinite dilution of MPCD-CC solutes as functions of the mass of solutes, for a dilute system with no direct interactions (N ed /L 3 box = 0.0013). Right: hydrodynamic radius a hyd deduced from Ds using the Stokes relation (Eq.( 9)).

hydrodynamic radius by changing the mass of the solute, which would modify the momentum exchange with the solvent particles during the collision step. We show in Fig. 2-left the self-diffusion coefficient at infinite dilution as a function of the mass M of the solute. Note that the results have been extrapolated to a simulation box of infinite size using eq. (15). We observe that the diffusion coefficient is an increasing function of the mass of the solute and reaches a plateau for large values of M/m f , more precisely for M larger than ≈ 2γ. This behavior was already observed by Ripoll et al 38 . Indeed, the heavier a solute is, the more it slows down solvent particles that participate in the collision step in the same cell. From these results, we extract the effective hydrodynamic radius from Stokes law (see Fig. 2-right). It appears then that the hydrodynamic radius can not exceed 0.35a 0 . We have also computed the self-diffusion at infinite dilution for a solute of mass M = 5m f in a MPCD solvent with γ = 10, keeping the same values of α and of δt c . We have obtained D • s = 0.02595 a 2 0 t -1 0 . For these parameters, the viscosity of the solvent is η = 8.7 m f a -1 0 .t -1 0 , so that the effective hydrodynamic radius in this case is a hyd = 0.23 a 0 . If we assume that the influence of the solute mass on the hydrodynamic radius is almost the same as the one described in Fig. 2, the case M = 5 m f , γ = 10 should be close to the case M = 3 m f , γ = 5. Indeed, we obtained a hyd ≈ 0.25 a 0 for M = 3 m f , γ = 5. Therefore, increasing the density of the MPCD fluid should not allow us to increase the hydrodynamic radius significantly.

Another route to determine the hydrodynamic radius of a solute is to induce a solvent flow around the solute from non-equilibrium simulations, and to compare the simulated flow to the analytical result, as described in Sec. II B. We have checked that the results were independent from the value of the solvent velocity in the small velocity regime. We have again studied the influence of the solute mass on its hydrodynamic radius. The mass of the particles varied from M = m f to M = 1000 m f . The Stokes flow field around a sphere with stick boundary conditions was well reproduced. The values of the hydrodynamic radius were fitted in order to minimize the difference between the computed and analytical velocity flow. We give in Fig. 3 the obtained hydrodynamics radii. In every case the root mean squared errors between computed and analytical velocities were between 1% and 2.5%. The results again show that the hydrodynamic radius is an increasing function of the mass which can be fitted by an exponential function, a hyd (M ) = 0.295 * (1 -exp(-M/3.3)). When the solute fixed at the center of the simulation box is heavy enough, the solvent particles in its cell are basically stopped because of the collision step, so that the hydrodynamic radius becomes independent from the solute mass at large mass. We also represent in Fig. 3 with a dashed line the exponential fit of the data obtained from equilibrium simulations (data of Fig. 2-right). Both methods, equilibrium and non-equilibrium simulations yield almost the same limiting value of the hydrodynamic radius at large solute mass: a hyd ≈ 0.3 a 0 . The fact that the hydrodynamic radius of solutes in MPCD-CC is in any case of the order of 0.3 a 0 is clearly a limitation of the method. For a given structural model of the solute (for instance the radius for a hard sphere model of solutes), the resolution of the MPCD grid relative to the size of the particles is imposed if one aims at keeping the structural radius close to the hydrodynamic one. Moreover, if the resolution of the grid is not thin enough, two or more solute particles can be located in the same MPCD collision cell at the same time. This raises the question of the validity of the method in such cases: What is the impact of the pathological description of short ranged hydrodynamic interactions when they emerge from momentum exchange including several solutes particles ? In the literature, the resolution of the grid is chosen so that the size of the collision cells roughly equals the diameter of the solute particles. This minimizes the probability that several solute particles participate in the collision step in the same cell. However, such choice leads to another kind of issue: The hydrodynamic radius is then different from the structural one. For a solute hard sphere diameter equal to 1.2 a 0 , the structural radius is indeed twice larger than the hydrodynamic one. This may lead to a decrease of hydrodynamic interactions, as we proceed to show.

IV. INFLUENCE OF THE SIZE OF THE COLLISION CELL ON THE DIFFUSION COEFFICIENT

We study in this part the influence of the resolution of the collision grid on the diffusion coefficient of solutes in collisional coupling with the solvent. We expect spurious effects on the dynamics of solutes when several solutes participate to the collision step in the same cell. These effects should be larger when the size of the cell increases and exceeds the minimal distance of approach between solutes, and/or when the concentration of solutes increases, and/or when attractions between solutes exist.

A. Solutes modelled as hard spheres

We have used the hard-sphere molecular dynamics algorithm 46 to compute the trajectories of solute particles between two collision steps in MPCD-CC. First, the influence of the volume fraction on the self-diffusion coefficient was studied for two different values of a 0 : a 0 = 0.4 a HS and a 0 = 3.0 a HS with a HS the hard sphere radius of the solute. For a 0 = 0.4 a HS the structural size is thus very close to the hydrodynamic one. The mass of the solute is M = 10 m f in every case. We present in Fig. 4 the computed self-diffusion coefficients divided by the value at infinite dilution (D • s = 4.22 10 -2 a 2 0 t -1 0 ) as functions of the volume fraction Φ HS = (4/3πN ed )(a HS /L box ) 3 . We have checked that the results obtained using a short-ranged Week-Chandler-Anderson interaction potential and a standard molecular dynamics algorithm instead of the hard-sphere algorithm coincide exactly.

As expected, the diffusion coefficient is in every case a decreasing function of the solute density (see Fig. 4). However, the influence of the size of the collision cell relative to that of the solute particles is striking. First, in the case a 0 = 3.0 a HS , the diffusion coefficients are considerably smaller than with a 0 = 0.4 a HS . Second, the decrease of the diffusion coefficient with the volume fraction is much more pronounced with a 0 = 3.0 a HS . This choice of the resolution of the grid is clearly not correct. We have thus more precisely investigated the influence of the size of the collision cell on the diffusion coefficient. We give in Fig. 5 the diffusion coefficients obtained at a given volume fraction as functions of the ratio between the size of the collision cell and the hard-sphere radius of solute a 0 /a HS . We focus on two volume fractions, φ HS = 0.05 and φ HS = 0.20. For small values of a 0 , i.e. for thin collision grids, the diffusion coefficient is almost independent from the value of a 0 (see the circles and the triangles symbols in Fig. 5). However, the diffusion coefficient of the solute starts to decrease dramatically for a 0 larger than 1.5 a HS . At the highest volume fraction, the value of D s /D • s drops by 60% when a 0 increases from 1.4 a HS to 3.0 a HS . This behavior is less pronounced in the dilute case, because the probability to have several solutes in the same collision cell is smaller: D/D • drops by 22% when a 0 increases from 1.4 a HS to 3.0 a HS . This strong decrease of the diffusion coefficient comes from the spurious coupling between solutes when they are in the same collision cell. Indeed, the maximum distance between two solutes in a collision cell is equal to √ 3a 0 ≈ 1.73 a 0 , so that we must have 2 a HS > 1.73 a 0 or a 0 /a HS < 1.16 to ensure that two hard-sphere solutes cannot be in the same collision cell. When a 0 is larger than this value, it may happen that two solutes participate in the collision step in the same cell. For a few values of a 0 /a HS in the dilute case, we have changed the ratio between the solute mass and the average density of the MPCD solvent in a cell (M = 5 m f with γ = 10). In this case, the infinite dilution diffusion coefficient is

D • s = 0.02595 a 2 0 t -1 0 .
We expected the problem of momentum exchange between solutes to be reduced in this case: If two solutes are in the same collision cell, their total mass stays of the order of that of the solvent also participating to the collision. Nevertheless, if the decrease of the diffusion coefficient when a 0 /a HS increases is smaller in these conditions, it is still clearly visible (square symbols in Fig. 5). To get more insight into these results, we compared the value of the diffusion coefficient with references obtained in a previous study from our group 30 . In this previous study, we compared the predictions of Brownian Dynamics (BD) simulations with hydrodynamic interactions at the Rotne-Prager level, BD without hydrodynamic interactions, and MPCD simulations with a central force coupling (CFC) between solutes and a MPCD fluid. Self-diffusion coefficients obtained from BD with hydrodynamics and from MPCD-CFC were in excellent agreement for neutral solutes interacting through a WCA potential. These results are independent from the size of the collision cell, and are represented in Fig. 5 as horizontal dashed lines for volume fractions equal to 0.05 and 0.20. We also report in the same figure the results obtained without hydrodynamic interactions as horizontal solid lines.

For the less concentrated system, φ HS = 0.05, the solid and dashed lines are close to each other because hydrodynamic interactions have a weak influence on the diffusion coefficient. The MPCD-CC results stand between the reference results with and without hydrodynamic interactions, which means that this method is unable to capture such subtle effect.

In the more concentrated case, hydrodynamic interactions have a non-negligible effect on the diffusion coefficient: D s /D • s is equal to 0.61 without HI and to 0.70 with HI. In this case, we observe a rather good agreement between the diffusion coefficient predicted by MPCD-CC and the reference for small a 0 values. D s /D • s computed by MPCD-CC is actually found to weakly increase at small a 0 values. This is due to the increase of hydrodynamic interactions because of the increase of a hyd in this range. Indeed, as shown in previous section, a hyd ≈ 0.3 a 0 . In the limit of small a 0 , the hydrodynamic radius is thus negligible compared to the structural radius a HS . At a 0 /a HS = 0.5 we have a hyd ≈ 0.15 a HS and at a 0 /a HS = 1.3, we have a hyd ≈ 0.39 a HS . This increase of a hyd leads to an increase of the hydrodynamic interactions. D s /D • s reaches a maximum at about a 0 /a HS = 1.3 but for larger a 0 values several solutes can interact via the collision step. It should be noted that the maximum value of D s /D • s computed by MPCD-CC is in excellent agreement with our reference. On the contrary, for larger values of a 0 , the MPCD-CC results tend towards those obtained without HI and becomes even smaller than them.

It is interesting to look at what happens for a 0 /a HS = 2, which is the most common choice in the literature when MPCD-CC is used [31][32][33] . For both volume fractions of 0.05 and 0.2, the value of the diffusion coefficient is below the one that was computed without hydrodynamic interactions. Does this mean that there are no hydrodynamic interactions ? No, since for this kind of parameters one can recover some features of hydrodynamic interactions, such as the scaling of polymer diffusion coefficients with polymer size. However, one should be careful when interpreting the results because an artificial coupling between solutes located in the same cell can happen in this range.

B. Solutes with electrostatic interactions

One may wonder whether the attraction between oppositely charged particles can change the impact of the size of the collision cell on the transport coefficients of solutes. To answer this question, we have computed the self-diffusion coefficients of a mixture of charged hard spheres of opposite charges, for several volume fractions between 0.033 and 0.185. The system corresponds to a 1-1 electrolyte solution described by the primitive model of electrolytes. The characteristic electrostatic length scale, the Bjerrum length l b = e 2 /(4πε 0 ε r k B T ) with e the elementary charge, ε 0 the permittivity of vacuum and ε r the relative permittivity of the solvent, is equal to 0.71 nm, corresponding to water at room temperature. We have compared the results for two different values of the cell size, a 0 /a HS = {0.4; 1.3}. For these grid resolutions, we have shown in the previous section that no spurious effect on the diffusion coefficient of neutral hard spheres was observed. As previously, we compare the self-diffusion coefficients to the reference values obtained in our previous paper i) with hydrodynamic interactions (MPCD with Central Force Coupling), and ii) without hydrodynamic interactions (Brownian Dynamics without HI). Results are shown in Fig. 6. They are close to those obtained with neutral hard-spheres: Diffusion coefficients increase slightly when the size of the collision cell increases, because the hydrodynamic radius increases, and this effect is more pronounced at large volume fractions. However, the computed values are not closer to the reference with hydrodynamic interactions than to the reference without hydrodynamic interactions. This contrasts with the case of neutral hard-spheres for φ = 0.2 and a 0 /a HS = 1.3 for which the computed diffusion coefficient was clearly closer to the reference with hydrodynamic interactions. The interplay between electrostatic and hydrodynamic effects may explain this discrepancy, since attractive coulombic interactions increase the impact of short range hydrodynamic interactions.

In the following, our goal is to compare MPCD-CC simulations with another reference that is known to describe accurately hydrodynamic interactions, with a ratio a hyd /a HS that is exactly the same as with MPCD-CC. In what follows, we have chosen to compute another transport coefficient of electrolyte solutions using MPCD-CC, the electrical conductivity, because a reliable semi-analytical theory of this quantity is available, which includes hydrodynamic interactions quantitatively.

V. ELECTRICAL CONDUCTIVITY OF AN IONIC SOLUTION

The measurement of the electrical conductivity is a widely used technique to analyze an electrolyte solution or a charged colloidal suspension known to be affected by hydrodynamic interactions. As ions of opposite charge move and drag the solvent in opposite directions under an electric field, hydrodynamic interactions between ions strongly reduce the conductivity compared to the infinite dilution value. Within the framework of the semianalytical theory described in Section II C, the hydrodynamic radius and the structural radius of ions can be treated as independent parameters. The structural radius a HS is a parameter of the interaction potential within the primitive model of electrolytes and is involved in the integral equations solved with the HNC closure to compute the pair distribution functions. Moreover, this radius is used as the lower bound of the integral which allows us to compute the electrostatic relaxation velocity correction δv elec i . The effective hydrodynamic radius of the charged solute which is related to the infinite dilution diffusion coefficient D • s through the Stokes law influences i) the electrical conductivity at infinite dilution, also called the ideal conductivity, ii) the electrostatic relaxation velocity correction δv elec i . When this theoretical framework is used to predict the electrical conductivity of real systems, a HS is usually obtained by mapping the equilibrium properties of the system to that of the primitive model of electrolytes, for example by fitting the osmotic coefficients 41,47 . The hydrodynamic size a hyd is usually deduced from the experimental asymptotic selfdiffusion coefficient at infinite dilution 48 . For most electrolytes, a HS and a hyd are close to each other 41,47 .

In our case, the comparison between MPCD-CC simulations and the theory will enable us to check whether both methods predict the same influence of hydrodynamic and electrostatic interactions on the electrical conductivity, for a given value of the hydrodynamic radius constrained by the MPCD-CC technique. In what follows, we study the same model of 1-1 electrolyte solution as in the previous section, with the size of the collision cell a 0 /a HS = 1.3. We stick to this value as it maximises the intensity of hydrodynamic interactions, as shown previously, without leading to spurious effects as two solutes cannot be in the same collision cell. MPCD-CC simulations were run for 2 10 5 t 0 . If the value of t 0 in real units is obtained by mapping the viscosity of water, it corresponds to a physical time of 23 ns. We have checked in every case that the radial distribution functions between ions obtained from MPCD-CC simulations were very close to those obtained within the HNC closure. The interest to use HNC total pair distribution functions in the theory instead of those computed from simulations is that they are obtained for larger distances and do converge to zero without statistical noise. We give in Fig. 7 the electrical conductivity as a function of the molar concentration of the electrolyte c = N ed N A /(2L 3 box ) with N A the Avogadro number. Error bars of the MPCD-CC results (open squares) are derived from the comparison of 10 independent simulations for each concentration. Results are compared to those obtained with the transport theory (filled circles). The electrical conductivity σ is divided by its value at infinite dilution σ • , which reads for solutions of charged solutes of opposite charge (1-1 electrolyte):

σ • = e 2 N ed D • s k B T L 3 box ( 16
)
with e the elementary charge.

The agreement between MPCD-CC and the theoretical predictions is excellent. For all systems but the more concentrated one, the values derived from the theory largely fall within the error bars of MPCD-CC results. For the system at 6.3 mol.L -1 , it is impossible to say whether the disagreement comes from a limitation of the theory or from the simulations. We also plot in Fig. 7 the results obtained with the theory without hydrodynamic interactions (filled squares and dashed lines): Only the correction to the velocity due to the electrostatic couplings, δv elec i , was taken into account in the theory. It is equivalent to suppress hydrodynamic interactions from the transport theory. We see in Fig. 7 that the values of the electrical conductivity are much higher in this case, and closer to the ideal conductivity (σ/σ • is close to 1). Moreover, we plot the theoretical results obtained by using a hyd = a HS instead of a hyd = 0.4 a HS (see the triangles in Fig. 7), i.e. in a case where hydrodynamic couplings are expected to be increased. The electrical conductivity is found to be strongly decreased compared to the case where a hyd = 0.4 a HS . In conclusion, these results show that : (i) MPCD-CC is able to capture hydrodynamic couplings corresponding to an hydrodynamic radius a hyd = 0.4 a HS , and (ii) the agreement between the theory with hydrodynamic couplings and the MPCD-CC requires that both treatments include hydrodynamic couplings with a quantitatively similar intensity.

VI. CONCLUSION

The use of mesoscopic simulation techniques to interpret dynamical properties of solutions usually requires special care, in particular when one is interested in a quantitative mapping between a real system and the simulated model. Traditional techniques, such as Brownian Dynamics, benefit from the fact that they are tightly related to the usual theories, and that the input parameters such as the infinite dilution diffusion coefficient or the hydrodynamic radius are widely used. Experimental studies have been improved to determine this class of parameters. More recent mesoscopic techniques, such as MPCD, rely on input parameters that are totally specific to the technique, and can be related to the properties of the fluid within a particular set of approximations. In order to understand the limitations of such methodologies, it can be useful to quantify commonly used dynamical quantities such as the hydrodynamic radius.

In this study, we find that the effective radius of a solute in collisional coupling with a MPCD solvent is about one third the size of the collision cell. To ensure that two hard-sphere solutes cannot be in the same collision cell, we should have in principle a 0 /a HS < 1.16 with a HS the hard-sphere radius of the solute. Moreover, we have shown that for neutral hard-spheres, hydrodynamic interactions between solutes are maximum at about a 0 /a HS = 1.3. In addition to these results, it should be noted that since the mean free path λ needs to stay close to 0.1 to keep the Schmidt number low 38 , the value of the kinematic viscosity ν is also constrained. Such constraints on the choice of MPCD parameters can be rephrased in term of time scales. When the solution is infinitely diluted, there are no hydrodynamic interactions. The only time scale of importance for transport is given by the infinite dilution diffusion coefficient, or by the product ηa hyd . The MPCD dynamics can be mapped onto a real time dynamics by mapping the MPCD diffusion coefficient to the value of the real diffusion coefficient of the solute. When solute concentration increases, a second important time scale emerges, the time scale at which hydrodynamic interaction are propagated, through the diffusion of momentum in the solution. This time scale is related to the kinematic viscosity ν. These two time scales are intimately related, as the diffusion coefficient depends on viscosity: At infinite dilution, we have D • s = k B T /(6πηa hyd ). The constraint on the hydrodynamic radius a hyd which depends on the size of the collision cell makes it impossible to have both D and ν match the values of a real system. For instance, let us consider a typical ion in water at ambient temperature, like a chloride or a sodium ion. Its diffusion coefficient is of the order of 10 -9 m 2 s -1 , and the dynamic viscosity of water is close to 10 -3 kg m -1 s -1 . The structural radius of an ion is about 0.15 nm. In our MPCD simulation, for typical values of the parameters (M = 10 m f , γ = 5, α = 130, λ = 0.1), the diffusion coefficient is equal to 4.22 10 -2 a 2 0 t -1 0 , and the kinematic viscosity is equal to 0.809 a 2 0 t -1 0 . The value of the MPCD time unit t 0 can be mapped into a real unit by choosing either the diffusion coefficient D, setting D M P CD = D real or η, setting ν M P CD = ν real . The two values of t 0 obtained by these mappings differ by one order of magnitude.

In the case of a polymer made of N monomers of infinite dilution diffusion coefficient D • s , within the Zimm theory, one can express the diffusion coefficient as a sum of a term depending on the diffusion coefficient, and a term depending directly on the viscosity and on an equilibrium property (the average inverse distance between monomers):

D = D • s N + k B T 6πη 1 r ij ( 17
)
The typical scaling appears when the first term becomes negligible, i. e. when the infinite dilution diffusion time scale is no longer relevant, and only hydrodynamic interactions affect transport. In such case, the main dynamic properties of the system should be qualitatively reproduced, and the MPCD time scale can be mapped onto the dynamics of the real system by setting ν M P CD = ν real . This is in agreement with the results obtained with polymers using the same simulation methodology. Although the hydrodynamic friction per monomer is very small, the correct scaling of diffusion properties on the number of monomers per polymer chain recovers the Zimm prediction, with 1 rij computed using the adequate theory of polymer equilibrium properties.

On the other hand, things might become more complicated when the scaling regime is not reached, for instance when finite size effects are important, typically for small values of N in the case of polymers. In the case of electrolytes, for most systems there is no regime in which hydrodynamic interactions totally dominate. In such case, although it is not possible to quantitatively map the dynamics of MPCD to that of a real system, the mesoscopic simulation can be used to assess the validity of a theory for which the parameters have been artificially chosen to be equal to those in MPCD. Such strategy is indeed the only possible way to explore the range of validity of theories for systems such as suspensions of nanoparticles or nanoporous media. Indeed, for such systems the interpretation of experiments often relies on concepts and parameters, such as the zeta potential, whose relationship with the underlying molecular system can be questionable 49,50 .

FIG. 3 :

 3 FIG. 3:Hydrodynamic radius a hyd of a solute in collisional coupling with a MPCD solvent determined from non-equilibrium simulations as a function of the solute mass M/m f . The solid line is an exponential fit of the data. The dashed line corresponds to an exponential fit of the data given in Fig.2-right.

FIG. 4 :

 4 FIG. 4: Diffusion coefficient of neutral hard spheres as function of their packing fraction φHS, using MPCD-CC simulations. The size of the collision cells is either a0 = 0.4 aHS, or a0 = 3.0 aHS.

FIG. 5 :

 5 FIG. 5:Diffusion coefficient of neutral hard spheres as function of the size of the collision cell a0 relative to the radius of the solute particle aHS, using MPCD-CC simulations. Two packing fractions are considered, φHS = 0.05 and a crowded case φHS = 0.20. The reference dashed lines correspond to MPCD with central force coupling, and the solid ones to Brownian Dynamics without hydrodynamic interactions.

FIG. 6 :

 6 FIG. 6:Self-diffusion coefficient of charged hard spheres (1-1 electrolyte) as function of their volume fraction φHS.The size of the collision cell is either a0 = 0.4 aHS, or a0 = 1.3 aHS. The reference dashed line corresponds to MPCD with central force coupling, and the solid one to Brownian Dynamics without hydrodynamic interactions.

FIG. 7 :

 7 FIG. 7: Electrical conductivity σ divided by the value at infinite dilution σ • , as a function of the electrolyte concentration. Open squares: MPCD-CC simulations with a0/aHS = 1.3; Filled circles: HNC-transport theory with the same parameters as MPCD-CC (D • s = 4.2210 -2 a 2 0 t -1 0 , aHS = a0/1.3); Filled squares: HNC-transport theory without hydrodynamic interactions; Filled triangles: HNC-transport theory with a hyd = aHS.

  . Self-diffusion coefficient of MPCD-CC solutes at infinite dilution as a function of the inverse box length. There are no direct interactions between solutes. The linear fit of the data writes Ds = 0.0422 -0.0382 1 L box .

			simulation result
		0.042	linear fit
		0.041
	]
	0 t 2 0 -1	0.04
	D s [a	0.039
		0.038
		0.037
		0 0.036	0.02 0.04 0.06 0.08	0.1	0.12
			1/L box [a 0 -1 ]
	FIG. 1:


Clément Ballabriga

Julien Forget

Laure Gonnord
email: laure.gonnord@ens-lyon.fr

Giuseppe Lipari

Jordy Ruiz

Vmcai

Static Analysis Of Binary Code With Memory Indirections Using Polyhedra

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In real time systems, checking that computations complete before their deadlines under all possible contexts is a crucial activity. Worst-Case Execution Time (WCET) analysis consists in computing an upper bound to the longest execution path in the code. It is usually performed on the binary code, because it needs information on the low-level instructions executed by the hardware processor.

In this paper, we propose a static analysis of binary code based on abstract interpretation using a polyhedra-based abstract domain. Our motivation is the need to enhance existing WCET analysis by improving the computation of upper bounds on the number of iterations in loops. However, our abstract domain has other potential applications (not developed in this paper), such as buffer-overflow analysis, unfeasible paths analysis or symbolic WCET computation [START_REF] Ballabriga | Symbolic wcet computation[END_REF].

Most analyses by abstract interpretation proposed in the literature are performed on source code. On the contrary, as it is usually the case for WCET analysis, we propose to analyze binary code. There are several important advantages in performing static analysis of binary code: 1) we analyze the code that actually runs on the machine, hence no need for additional assumptions on how the compiler works; 2) in presence of undefined behaviors (of source code), the analysis is more accurate; 3) we can perform the analysis even without access to the source code.

The main problem is that, in higher-level representations, the variables, addresses and values are well identified. In binary code, the notion of program variable is lost, so we can only analyze processor registers and memory locations. We propose to identify the subset of registers and memory locations to be represented in the abstract state as the analysis progresses. This representation enables us to design a relational analysis on binary code, which is the main contribution of the paper.

Motivating example

As a motivating example, we present a snippet of C code, inspired from packet processing network drivers in Figure 1 3 . We remind however that our methodology addresses (disassembled) binary code.

The send_request function sends a request in some application-layer protocol that runs over UDP/IP. Lines 12-13 build a packet composed of a variablelength IP header, a fixed-length UDP header, and a variable-length UDP payload (some operations on IP or UDP fields have been omitted). Note that the starting address of the UDP header depends on the size of the IP header (h1->hdr_len). At line 17, we call the function responsible for putting the useful data (payload) into the packet. At line 18, the packet is sent using the send_packet function, which belongs to the lower-level network layer API. This function does not take the packet size as parameter, since it can be deduced from the header: in lines 2-3, the function parses the packet to obtain the UDP payload size, and the UDP checksum is computed by iterating over the payload.

To automatically compute a bound on the number of iterations of the loop at line 4, the analysis has to discover that udp_l equals udp_size (due to line 16). This can be done with an appropriate use of a relational abstract domain. However, very few of the existing analyses running on binary code use a relational domain, and to the best of our knowledge, none support relations between addresses that are not know statically (udp_l, udp_size). Let us emphasize that such a use of pointers and memory buffers is typical of many embedded systems: for instance in network packet processing, but also in many device drivers.

Contribution

The contributions of the paper are:

-A new relational abstract domain POLYMAP, which consists of a polyhedron and two mappings that track the correspondence between data locations (registers or memory) and polyhedra variables; -An abstract interpretation procedure, which computes abstract states of POLYMAP for a small assembly language, and which we prove to be sound; -An experimental evaluation of our prototype called Polymalys. It implements the previous procedure and computes upper bounds to loop iterations. We compare Polymalys with other existing tools on a set of classic benchmarks.

Language definition

In this section, we define the analyzed language, called MEMP, a simplified assembly language where we focus on memory indirection operators.

Syntax

In order to simplify the presentation, we make the following assumptions: all data locations have the same size, memory accesses are aligned to the word size, there are no integer overflows, and function calls are inlined (these limitations could be lifted using for instance [START_REF] Bygde | Fully bounded polyhedral analysis of integers with wrapping[END_REF][START_REF] Sharir | Two approaches to interprocedural data flow analysis[END_REF]). We also reduce the set of instructions to a minimum (Polymalys actually supports the ARM A32 instruction set). The syntax of MEMP is defined in Figure 2. A program is a sequence of labeled instructions. Instructions operate on registers, labels or constants. Concerning memory instructions, if r contains value c, then * (r) denotes the content at address r (below, we overload the notation and also denote * (c) for this content).

OP c denotes the concrete semantics of operation OP. RAND emulates undefined registers, to represent e.g. function parameters. Other instructions are directly commented in the figure (on the left of each instruction). :

:= r1 ← OP c (r2, r3) | OP r1 r2 r3 r ← c | SET r c Emulate undefined r | RAND r r1 ← * (r2) | LOAD r1 r2 * (r1) ← r2 | STORE r1 r2 Branch to l if r = 0 | BR r l Halt | END Fig. 2. Syntax of MEMP

Formal semantics

The small-steps semantics of MEMP is defined below. The semantics of data and arithmetic/logic operations is defined in Figure 3 by function i -→, which operates in a context (R, * ) consisting of two mappings where:

-R : R → Z is the registers content, which maps registers to their values. We assume that it is initially empty; - * : Z → Z is the memory content, which maps memory addresses to their values. We assume that it is also initially empty. Note that integer wrapping could be used to restrain addresses to be in N instead of Z [START_REF] Bygde | Fully bounded polyhedral analysis of integers with wrapping[END_REF].

For a given mapping m, we denote m[x : y] the mapping m such that m (x) = y and, for every register x = x, m (x ) = m(x ). In other words, m[x : y] denotes a single mapping substitution (or mapping addition if x was previously unmapped). We also denote m\(x 1 : x 2 ) the mapping such that the association

x 1 : x 2 is removed from m.
The semantics of control flow operations is defined in Figure 4, by the function c →, which adds a program counter pc to the previous context. We use e → to denote the last transition of the program.

Abstract domain

The abstract domain we propose is based on the polyhedral abstract domain [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], to which we add information to track relations between polyhedra variables and registers or memory addresses.

(SET r c, R, * ) i -→ (R[r : c], * ) c = random() (RAND r, R, * ) i -→ (R[r : c], * ) R(r2) = c2 R(r3) = c3 c1 = OP c (c2, c3) (OP r1 r2 r3, R, * ) i -→ (R[r1 : c1], * ) R(r2) = c2 * (c2) = c1 (LOAD r1 r2, R, * ) i -→ (R[r1 : c1], * ) R(r1) = c1 R(r2) = c2 (STORE r1 r2, R, * ) i -→ (R, * [c1 : c2])
Fig. 3. Semantics of data and arithmetic operations.

P [pc] = BR r l R(r) = 0 P (pc, R, * ) c → (pc + 1, R, * ) P [pc] = BR r l R(r) = 0 P (pc, R, * ) c → (l,

Polyhedra

A polyhedron p denotes a set of points in a Z vector space bounded by linear constraints (equalities or inequalities). More formally, let |S| denote the cardinality of set S. Let C n denote the set of linear constraints in Z n on the set of variables V n , where |V n | = n. Then c 1 , c 2 , ..., c m denotes the polyhedron p consisting of all the vectors in Z n that satisfy constraints c 1 , c 2 , . . ., c m , where c i ∈ C n for 1 ≤ i ≤ m (n and m are unrelated). We denote dim(p) = n the dimension of p. In the rest of the paper, the term variable implicitly refers to polyhedron variables. We denote:

-P the set of polyhedra; s ∈ p when s (with s ∈ Z dim(p) ) satisfies the constraints of polyhedron p; p p iff ∀s ∈ p, s ∈ p ; p = p p the convex hull of p and p ; p = p p the union of the constraints of p and p ; vars(p) the set of variables of p, where |vars(p)| = dim(p) by definition; proj (p, x 1 . . . x k ) the projection of p on space x 1 . . . x k , with k < |dim(p)|; p[x i /x j ] the substitution of variable x j by x i in p.

Abstract States

In polyhedral analysis of source code, variables of the polyhedra are related to variables of the source code. In our case, polyhedra variables are related to registers and memory contents. We use the term data location to refer indistinctly to registers or memory addresses. Let V denote the set of polyhedra variables.

The set of abstract states POLYMAP is defined as

A = P × (R → V) × (V → V).
An abstract state a ∈ A, with a = (p, R , * ), consists of a polyhedron p, a register mapping R and an address mapping * . We have R (r) = v iff variable v represents the value of register r in p. We have * (x 1 ) = x 2 iff variable x 2 represents the value at the memory address represented by variable x 1 . We denote vars R (p) the codomain of R (i.e. register content variables), vars A (p) the domain of * (i.e. address variables) and vars C (p) the codomain of * (i.e. address content variables). Sets vars R (p), vars A (p) and vars C (p) are disjoint and are all subsets of vars(p).

Example 1. In the following abstract state, register r 0 contains value 2, and address 2 contains value 1:

({x 1 = 2, x 2 = x 1 , x 3 = 1}, {r 0 : x 1 }, {x 2 : x 3 })
The usual operators on the abstract domain (inclusion, join and widening), and its least and greatest elements are presented in Section 4.4.

Aliasing

In a general sense, aliasing occurs in a program when a data location can be accessed through several symbolic names. As we will see in Section 4, aliases play an important role in our analysis. In fact, we introduce mechanisms that prevent their occurrence in the abstract state (see Section 4.2), so as to simplify the analysis. We define the aliasing relation between two variables x 1 and x 2 of a polyhedron p as follows:

-Cannot alias: whenever

x 1 = x 2 ∩ p = ∅; -May alias: whenever x 1 = x 2 ∩ p = ∅; -Must alias, denoted x 1 ≡ x 2 : whenever p x 1 = x 2 .
The aliasing relation between a register r and a variable x is defined by the aliasing relation between R (r) and x. Similarly, the aliasing relation between two registers r 1 , r 2 is defined by the aliasing relation between R (r 1 ) and R (r 2 ).

To avoid ambiguities with notations on constraints, let same(x 1 , x 2 ) denote the fact that x 1 and x 2 are the same polyhedron variables (not just equivalent variables). There is no need to check register aliases, because a single register cannot be mapped to two different variables (R is a function). The absence of aliases can thus be stated as follows.

Definition 1. Let s = (p, R , * ) be an abstract state. We say that s is alias free iff:

∀x 1 , x 2 ∈ vars A (p), x 1 ≡ x 2 ⇒ same(x 1 , x 2 )

Computing abstract states

Our analysis follows the abstract interpretation framework proposed in [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], adapted to our setting with non-local control-flow, following the technique proposed in Astrée [START_REF] Kästner | Astrée: Proving the absence of runtime errors[END_REF] and MOPSA [START_REF] Miné | Design of a Modular Platform for Static Analysis[END_REF]. An important singularity of our analysis is that polyhedral variables are progressively created or removed during the analysis. Whenever a new polyhedron variable is introduced, we assume it is a fresh variable that has never been used at any other point during the analysis.

Interpretation algorithm

We use (p , [r i :

x i ], [x j : x k ])(•) as a shorthand for λ(p, R , * ).(p p , R [r i : x i ], * [x j : x k ]
), and denote -when a state component remains unchanged. Procedures to compute the join ( ), widening ( ) and antialias of abstract states, and the transfer function (I ) of instruction I are detailed in the remainder of this section. The complete interpretation procedure is described in Algorithm 1. It applies to a program P of MEMP. During the interpretation, we keep a subset L of labels of interest. Abstract values are stored in a map M from labels to abstract values. We assume that loop header labels L W of P have previously been identified using an existing analysis (e.g. Tarjan's algorithm [START_REF] Tarjan | Depth-first search and linear graph algorithms[END_REF]). Figure 5 reports a running example of this analysis, that will be used throughout the rest of the section.

Algorithm 1 Interpret(P)

1: procedure update( , a, L)

Auxiliary procedure 2:

a ← antialias(a)

3: if ∈ LW then Check if l is a loop header 4: new ← M [l ] (M [l ] a 5: else 6: new ← M [ ] a 7: end if 8: if new M [ ]
5 p1 = x1 = 4, x2 = 5 R 1 = {r0 : x0, r1 : x1, r2 : x2, r7 : x7} 6 p2 = p1 x3 = x0 + x1 R 2 = R 1 [r3 : x3] 7 p3 = p2 x4 = x3, x5 = x1 R 2 * 1 = {x4 : x5} 8 p4 = p3 x8 = x7 -x1 R 3 = R 2 [r5 : x8] * 1 10 (from 9) p5 = p4 x9 = x2 R 3 * 2 = {x4 : x9} 10 (from 8) p6 = p4 x8 = 0 R 3 * 1 unify(10, 10 ) p7 = p6[x9/x5] R 3 * 3 = {x4 : x9} 10 10 p8 = p2 x4 = x3, x8 = x7 -x1, R 3 * 3 x1 ≤ x9 ≤ x2 11 p8 x10 = x9 R 3 [r6 : x10] * 3
Fig. 5. Running example of analysis

Anti-aliasing

Whenever updating an abstract state, we immediately remove aliases (line 2), because the absence of aliases significantly simplifies the analysis in places where we need to check the equivalence of two variables (LOAD, STORE, and ). In practice, aliases are introduced when encountering a conditional branching (see Section 4.3). We remove an alias using procedure antialias, which relies on the procedure Merge defined below. It is based on the following observation: if two addresses are equal, then the values stored at these addresses must be equal too. Let x 1 , x 2 be two variables of vars A (p) such that:

¬same(x 1 , x 2 ) ∧ x 1 ≡ x 2 . Merge((p, R , * ), x 1 , x 2 ) = (p , R , * ) with p = (p[x 1 /x 2 ])[ * (x 1 )/ * (x 2 )] and * = * \ (x 2 : * (x 2 ))
Function antialias : A → A applies Merge for each pair of distinct equivalent address variables of an abstract state.

Example 2. In state a below, address x 2 is an alias on address x 1 . Thus, x 4 must be equal to x 3 , so Merge(a, x 1 , x 2 ) replaces x 2 by x 1 and x 4 by x 3 . In the result, x 3 is constrained by the original constraints of x 3 and x 4 , and the memory mapping x 2 : x 4 is discarded.

a = ( x 1 = x 2 , x 3 ≥ 4, x 4 ≤ 5 , -, * = {x 1 : x 3 , x 2 : x 4 }) Merge(a, x 1 , x 2 ) = ( 4 ≤ x 3 ≤ 5 , -, * = {x 1 : x 3 })

Transfer functions

We now define the constraints generated for the analysis of each instruction of our language. We denote (I ) : A → A the transfer function of instruction I.

Binary operation If the relation r 1 = OP c (r 2 , r 3 ) is linear, we map the target register to a new variable, subject to the corresponding linear constraint in the polyhedron. The memory mapping is unchanged. Otherwise, the target register is mapped to a new unconstrained variable.

(OP r 1 r 2 r 3 ) = ( x = OP c (R (r 2 ), R (r 3 )) , [r 1 : x], -)(•) if linear(OP c ) (-, [r 1 : x], -)(•) otherwise
Example 3. In Figure 5, at label 6 (i.e. the label immediately following the ADD operation) we introduce the constraint x 3 = x 0 + x 1 and the register mapping

R 1 (r 3 ) = x 3 .
Set The impact of the immediate load instruction is straightforward:

(SET r 1 c) = ( x = c , [r 1 : x], -)(•)
Rand The random instruction maps a register to an unconstrained variable:

(RAND r 1 ) = (-, [r 1 : x], -)(•)
Load If the input state contains a memory address variable that is equivalent to the load address (note that for alias free states, if such a variable exists, it is unique), then in the output state the value of the destination register is the value of the memory value mapped to this address. Otherwise, the value of the destination register is undefined:

(LOAD r 1 r 2 ) = ( x = * (a) , [r 1 : x], -)(•) if a ≡ r 2 (-, [r 1 : x], -)(•) otherwise
Example 4. In Figure 5, at label 10 we have x 4 ≡ r 3 and * (x 4 ) = x 9 , so at label 11 we introduce the constraint x 10 = x 9 and the mapping R 3 [r 6 ] = x 10 .

Store Again, we need to consider the impact of aliases. If there exists an address variable equivalent to the target register, then there already exists a memory mapping for this address.

Replace(a) = ( x = R (r 2 ) , -, [a : x])(•) Create = ( x i = R (r 1 ), x j = R (r 2 ) , -, [x i : x j ])(•) May = {a∈A|a may-alias r1}
λs.(Replace(a)(s) s)

Example 5. In Figure 5, at label 7, we create a new memory mapping * 1 (x 4 ) = x 5 and we introduce the constraints x 4 = x 3 , x 5 = x 1 .

Example 6. In Figure 5, at label 10, when coming from label 9, we replace a previous mapping, x 4 is mapped to x 9 (instead of x 5 previously), and we introduce the constraint x 9 = x 2 .

Branching In Algorithm 1, when branching to a target label ( ) the branching condition holds (r = 0). We add no constraint for the otherwise case because it cannot be encoded using a linear relation.

Example 7. In Figure 5, at label 10, when coming from label 6, we add the constraint x 8 = 0.

Abstract domain operators, least and greatest elements

Our analysis introduces new variables and removes old ones as it progresses.

There is no predefined correspondence between variables and data locations, because the set of data locations used by the program is unknown a priori. As a consequence, it may happen that two abstract states use different variables to designate the same data location. This implies that to compare two states we first need to check whether some variables of the two states actually correspond to the same data location. This verification relies on a unification procedure, presented below. Unification is used for inclusion testing, and also in the join and widening operators.

Unification Unification checks for the equivalence of two variables in two polyhedra, p 1 and p 2 . Intuitively, we try to express each variable as a linear expression of a well-chosen set of variables to conveniently check their equivalence.

Let V c = vars(p 1 ) ∩ vars(p 2 ) and p = proj (p 1 , V c ) proj (p 2 , V c ). We denote npiv (p ) the set of non-pivot variables discovered by Gauss-Jordan elimination performed on the system of equality constraints of p (we exclude inequalities). Then, npiv (p ) is such that, in p :

no variable in npiv (p ) is equivalent to a linear expression of other variables of npiv (p ); each variable in vars(p )\npiv (p ) is equivalent to a linear expression of variables from npiv (p ).

Let linexpr (x, p 1 , npiv (p )) denote the linear expression representation of variable x ∈ vars(p 1 ) in terms of variables in npiv (p ), represented as the vector of the linear expression coefficients. Let C be the constraint system of proj (p 1 , x∪npiv (p )). If C contains an equality constraint involving x, then computing linexpr (x, p 1 , npiv (p )) is straightforward. Otherwise, the empty vector is returned. If several (non-equivalent) equality constraints appear, we arbitrarily pick one. Note that, even though our unification can miss equivalent variables, this does not jeopardize the soundness of the analysis (see Section 5.3 and in particular Lemma 3).

Algorithm 2 describes our unification procedure. We directly modify the second state to unify it with the first one. First, we compute set of non-pivot variables (line 4). Then, we check for the equivalence of address variables according to their linear expression representation, and we perform variable substitutions in p 2 , R 2 and * 2 in case of equivalence (line 8). Register unification is simpler, we just replace the bindings in R 2 by those of R 1 (line 12). Example 8. In Figure 5, when computing unify(10, 10 ), s 1 corresponds to the state of 10 and s 2 to the state of 10 . A possible set of non-pivot variables is {x 0 , x 7 }. In s 1 (and in s 2 ), we have x 4 -x 0 + 0 • x 7 -4 = 0, so linexpr (x 4 ) = [1; -1; 0; -4] (corresponding, respectively, to the coefficients of x 4 , x 0 , x 7 , and the constant). Since * 2 (x 4 ) = x 9 (in s 1 ) and * 1 (x 4 ) = x 5 (in s 2 ), we replace x 5 by x 9 in s 2 .

Algorithm 2 unify((p

1 , R 1 , * 1 ), (p 2 , R 2 , * 2 )) 1: (p 2 , R 2 , * 2 ) ← (p2,
Inclusion Let us now define formally the partially ordered set (A, ). Given two functions f and g, we denote f ⊆ g when Dom(f ) ⊆ Dom(g) and ∀x ∈ Dom(f ) : f (x) = g(x). Introducing new mappings in R or * (i.e. enlarging their domains) actually removes feasible concrete states, thus we define abstract states inclusion as follows (see lemma 4 for more details):

Definition 2. Let a 1 = (p 1 , R 1 , * 1 ) and a 2 = (p 2 , R 2 , * 2 )
. The ordering operator is defined as follows:

a 1 a 2 ⇔p 1 p 2 ∧ R 2 ⊆ R 1 ∧ * 2 ⊆ * 1 with (p 1 , R 1 , * ) = unify(a 2 , a 1 )
There exists several equivalent representations of the greatest and least elements of (A, ). We define them as follows: Definition 3. The greatest element of (A, ) is denoted , with = ( , ∅, ∅). Definition 4. The least element of (A, ) is denoted ⊥ and defined as ⊥ = (p ⊥ , R ⊥ , * ⊥ ), where p ⊥ is the empty polyhedron and R ⊥ , * ⊥ are such that every data location is mapped to a variable.

Join Algorithm 3 describes our join procedure. It unifies the input states (line 1), then computes the convex hull on the unified states (line 2). Then, if a memory location or register is bound in one input state and unbound in the other, it is unbound in the result state. Example 9. In Figure 5, when computing 10 10 , we obtain identical register and memory mappings for 10 and unify [START_REF] Bygde | Fully bounded polyhedral analysis of integers with wrapping[END_REF][START_REF] Bygde | Fully bounded polyhedral analysis of integers with wrapping[END_REF]. The convex hull p 5 p 7 groups the constraints on x 9 (x 1 ≤ x 9 ≤ x 2 ) and lifts those on x 8 .

Algorithm 3 (p 1 , R 1 , * 1 ) (p 2 , R 2 , * 2 ) 1: (p 2 , R 2 , * 2 ) = unify((p1, R 1 , * 1 ), (p2, R 2 , * 2 )) 2: p ← p1 p 2 3: R ← ∅; * ← ∅ 4: for all r ∈ Dom(R 1 ) do 5: if R 1 (r) = R 2 (r) then R (r) ← R 1 (
Widening Due to the presence of loops, the widening operator is used to ensure that our analysis reaches a fixpoint. is defined just like , except that we use a polyhedra widening operator in place of .

Loop bounds

To compute loop bounds, for each loop header label we create a "virtual" register r , to count the number of iterations of . We instrument the program so that the register r is set to 0 when entering loop , and incremented at each iteration of (which is fairly classic, see e.g. [START_REF] Gulwani | Speed: precise and efficient static estimation of program computational complexity[END_REF]). Finally, let P a program of MEMP and M = interpret(P ). Let e be the label of instruction

END in P . Let (p f , R f , * f ) = M [ e ].
Then the loop bound for a loop header is computed as max(p f , R f [r ]) (where max(p, x) denotes the greatest value of variable x satisfying the constraints of p).

Soundness

In this section, we prove the soundness of our analysis. We first establish a set of important lemmas on our abstract domain operators, and then prove soundness with respect to the concretization function.

Join

Operator is not commutative. We establish that it does however compute an upper bound of its operands, with respect to our inclusion definition (Lemma 1). The proof is based on two auxiliary properties on mapping inclusions:

Property 1. Let a 1 = (p 1 , R 1 , * 1 ), a 2 ∈ A, a 3 = (p 3 , R 3 , * 3 ) = a 1 a 2 . We have: (p 1 p 3 ) ∧ (R 3 ⊆ R 1 ) ∧ ( * 3 ⊆ * 1 )
Proof. Considering Algorithm 3: (p 1 p 3 ) follows from line 2, (R 3 ⊆ R 1 ) from line 5, and ( * 3 ⊆ * 1 ) from line 8.

Property 2. Let a 1 , a 2 , a 1 ∈ A, with a 1 = (p 1 , R 1 , * 1 ) = unify(a 2 , a 1 ). Then: (R 2 ⊆ R 1 ) ∧ ( * 2 ⊆ * 1 ) ⇒ (R 2 ⊆ R 1 ) ∧ ( * 2 ⊆ * 1 ) Proof. Obvious from Algorithm 2. Lemma 1. Let a 1 , a 2 ∈ A. We have: (a 1 a 1 a 2 ) ∧ (a 2 a 1 a 2 ).
Proof. Polyhedron inclusion follows from the polyhedra join operator. We must also prove the inclusion of register and memory mappings (after unification).

Case for a 1 follows from Properties 1 and 2. Concerning the case for a 2 , let a 3 = a 1 a 2 . When computing a 3 , a variable v of a 2 falls into one of three categories: 1) v is also in vars(p 1 ), it remains in a 3 ; 2) v is equivalent to a variable v 1 of vars(p 1 ), it is replaced by v 1 in a 3 (Algorithm 2, line 8); 3) otherwise, it is removed (Algorithm 3). Then, let a 2 = unify(a 3 , a 2 ). When computing a 2 , variables that fell in category 2 at the previous step (when computing a 3 ) will be replaced by their equivalent in a 3 , because they fall again in category 2. Thus we obtain R 3 ⊆ R 2 , * 3 ⊆ * 2 , which concludes the proof.

Widening

Lemma 2 establishes that operator is indeed a widening operator.

Property 3. Let a 1 , a 2 ∈ A. We have: (a 1 a 2 ) (a 1 a 2 ).
Proof. The property holds because and use the same unification procedure, and because we assume that is a valid polyhedra widening operator.

Property 4. Let a 1 = (p 1 , R 1 , * 1 ), a 2 ∈ A, a 3 = (p 3 , R 3 , * 3 ) = a 1 a 2 . We have: (p 1 p 3 ) ∧ (R 3 ⊆ R 1 ) ∧ ( * 3 ⊆ * 1 )
Proof. Same as for Property 1.

Property 5. Let (b n ) n∈N be a non decreasing infinite sequence in A. Then, the sequence a 0 = b 0 and a n+1 = a n b n+1 converges in a finite number of steps.

Proof. Thanks to Property 4, and considering that there is a finite quantity of data locations, there exists N ∈ N such that for all i > N , R i+1 = R i and * i+1 = * i . Thus, a i+1 = (p i q i+1 , R i , * i ), where q i+1 is the polyhedron of b i+1 and p i that of a i .

Assuming that is a valid polyhedra widening operator, there exists m > N such that p m+1 = p m . Since m > N we also have R m+1 = R m and * m+1 = * m , which concludes the proof. Lemma 2. Operator is a widening operator.

Proof. Follows from Properties 3 and 5.

Concrete and abstract states

Let C = ((R → Z) × (Z → Z)) denote the set of concrete states (pairs of registers contents and memory contents). Data locations are mapped to values in a concrete state, while they are mapped to polyhedra variables in the abstract state. The concretization function γ relates data location values to data location variables as follows:

Definition 5. Let a = (p, R , * ) be an abstract state. The concretization function γ is defined as follows:

γ : A -→ P(C) (p, R , * ) -→ ( * , R) | ∃f : Dom( * ) → Dom( * ),   r∈Dom(R ) R (r) = R(r) x∈Dom( * ) x = f (x), * (x) = * (f (x))   p
More intuitively, we build a polyhedron p with the following constraints: 1) register values of the concrete state (R(r)) must be equal to the corresponding variable in the abstract state (R (r)); 2) we try to find a function f that maps address variables to addresses (x = f (x)), then the content of each address variables ( * (x)) must be equal to the memory value ( * (f (x))). If p p then the concrete state satisfies the constraints of p and belongs to the concretization.

Example 10. a = ({1 ≤ x 1 ≤ 2, x 2 = x 1 , x 3 = 1}, {r 0 : x 1 }, {x 2 : x 3 }) γ(a) = {({r 0 = 1}, { * (1) = 1}}), (f (x 2 ) = 1) ({r 0 = 2}, { * (2) = 1})} (f (x 2 ) = 2) Let c - → * denote the transitive closure of c →.
The soundness of our abstract interpretation is established as follows:

Theorem 1. Let P be a MEMP program. Let M = Interpret(P ). Then, for any concrete state s init :

(P (l 1 , s init ) c - → * ( , s)) =⇒ (s ∈ γ(M [ ]))
Proof. The proof of soundness follows from the structure of Algorithm 1, and from the following lemmas, which establish the soundness of each operator used in the algorithm. Lemma 3. Let a 1 , a 2 ∈ A. We have: γ(a 1 ) = γ(unify(a 2 , a 1 )).

Proof. Let a 1 = unify(a 2 , a 1 ). Since we assume that a 1 and a 2 are alias free (recall Section 4.2), any two non-equivalent variables in a 1 are also replaced by non-equivalent variables in a 1 (or unchanged). Thus a 1 is a simple renaming of a 1 , and so a 1 and a 1 have the same concretization.

Lemma 4. Let a 1 , a 2 ∈ A. We have: (a 1 a 2 ) ⇒ γ(a 1 ) ⊆ γ(a 2 )
Proof. Let s ∈ γ(a 1 ). Let a 1 = (p 1 , R 1 , * 1 ) = unify(a 2 , a 1 ). From Lemma 3, s ∈ γ(a 1 ), thus there exists a function f for s satisfying the property of Definition 5 with a = a 1 . Now, assume that p 1 p 2 ∧ R 2 ⊆ R 1 ∧ * 2 ⊆ * 1 (i.e. a 1 a 2 ). Then there exists a function f for s that satisfies Definition 5, with a = a 2 : just take f such that it is the restriction of f to Dom( * 2 ). So s ∈ γ(a 2 ).

Several other abstract domains other that polyhedra, capable of representing linear constraints between variables, have been proposed, such as for instance [START_REF] Karr | Affine relationships among variables of a program[END_REF][START_REF] Venet | The gauge domain: Scalable analysis of linear inequality invariants[END_REF][START_REF] Miné | Field-sensitive value analysis of embedded C programs with union types and pointer arithmetics[END_REF]. Choosing the most appropriate domain boils down to a trade-off between the execution time and the precision of the analysis. In our work we chose the polyhedra domain and thus favored precision. However, we think that it would be simple to adapt our work to another domain (e.g. to reduce analysis time), because our computation of memory and register mappings does not depend on how constraints between variables are represented and computed.

Several works address static analysis of binary code [START_REF] Balakrishnan | Analyzing memory accesses in x86 executables[END_REF][START_REF] Djoudi | Binsec: Binary code analysis with low-level regions[END_REF][START_REF] Reps | Improved memory-access analysis for x86 executables[END_REF][START_REF] Sepp | Precise static analysis of binaries by extracting relational information[END_REF][START_REF] Bardin | Refinement-based CFG reconstruction from unstructured programs[END_REF], however they do not consider the problem of identifying memory locations of interest. In contrast, we identify these locations during the analyses.

An important problem when dealing with binary code analysis is to figure out the set of interesting data locations used by the program. This is related to pointer analysis (the so-called aliasing problem), and has been extensively studied [START_REF] Hind | Pointer analysis: Haven't we solved this problem yet?[END_REF][START_REF] Hardekopf | The ant and the grasshopper: fast and accurate pointer analysis for millions of lines of code[END_REF]. While the majority of pointer analyses have been proposed in the context of compiler optimizations, a certain number of ideas can be borrowed and applied to binary code analysis.

In this paper, our approach is applied to static loop bound estimation, in the context of WCET analysis, so we compare our results with other loop bound estimation tools. The oRange tool [START_REF] Bonenfant | oRange: A tool for static loop bound analysis[END_REF] is based on an abstract interpretation method defined in [START_REF] Ammarguellat | Automatic recognition of induction variables and recurrence relations by abstract interpretation[END_REF]. It provides a very fast estimation of loop bounds, but it is restricted to C source code. SWEET [START_REF] Ermedahl | Loop bound analysis based on a combination of program slicing, abstract interpretation, and invariant analysis[END_REF] features a loop bound estimator, which works on an intermediary representation (ALF format). The approach is based on slicing and abstract interpretation and it generally provides very tight loop bounds even in complex cases, but the running time of the analysis seems to depend on the loop bound values, and in our experience for large loop bounds the analysis did not terminate.

KTA [START_REF] Broman | A brief overview of the KTA WCET tool[END_REF] is a static WCET analysis tool based on abstract interpretation and path exploration of binary code. As its purpose is to compute a WCET, it does not directly provide information on loop bounds and we could not find documentation on the method used to compute these bounds. Thus, KAT was not included in our benchmarks. Furthermore, the analysis time seems to depend on the loop bound values.

Compared to these existing works, our approach combines the polyhedral domain with binary code analysis, taking into account memory accesses and supporting analysis of relations between unknown memory addresses; moreover our method is proved to be sound and to always terminate.

Experimental results

Our methodology is implemented in a prototype called Polymalys. Our experiments consist of two parts. First, we validate our approach by comparing Polymalys with other existing loop bound analysis tools on classic benchmarks. Then, we provide detailed examples of programs for which Polymalys successfully estimates loops bounds, while the other tools fail to do so.

Implementation

Polymalys is implemented as a plugin of OTAWA (version 2.0), an open source WCET computation tool [START_REF] Ballabriga | OTAWA: An open toolbox for adaptive WCET analysis[END_REF]. Polymalys relies on OTAWA for control-flow analysis and manipulation, and on PPL [START_REF] Bagnara | The Parma polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] for polyhedra operations. Polymalys implements several optimizations to reduce the number of variables and constraints of an abstract state (p, R , * ), most notably:

-Unmapped variables: any variable that is not in R or in * can be safely removed from the polyhedron by performing a projection on the remaining (used) variables; -Dead registers: we remove dead register variables by perform a preliminary liveness analysis, using classic data-flow analysis methods [START_REF] Aho | Compilers: Principles, Techniques, and Tools[END_REF]; -Out-of-scope variables: whenever modifying the stack pointer register (SP ), assuming that the stack grows downwards, for each pair of variables (x i , x j ) such that * (x i ) = x j , if p x i < R [SP ] then x i and x j can be removed.

Benchmarks

The analyses have been executed on a PC with an Intel core i5 3470 at 3.2 Ghz, with 8 GB of RAM. Every benchmark has been compiled with ARM crosstool-NG 1.20.0 (gcc version 4.9.1), using the -O1 optimization level. First, we report the results of our experiments on the Mälardalen benchmarks [START_REF] Gustafsson | The Mälardalen WCET benchmarks: Past, present and future[END_REF] and on PolyBench [START_REF] Pouchet | Polybench: The polyhedral benchmark suite[END_REF] in Table 1. The benchmarks gemver, covariance, correlation, nussinov and floyd-warshall are from PolyBench, while the others are from Mälardalen. We exclude benchmarks that are not supported by OTAWA, mainly due to floating point operations or indirect branching (e.g. switch). We compare Polymalys with SWEET [START_REF] Lisper | SWEET-a tool for WCET flow analysis[END_REF], PAGAI [START_REF] Henry | Pagai: A path sensitive static analyser[END_REF] and oRange [START_REF] Bonenfant | oRange: A tool for static loop bound analysis[END_REF]. For each benchmark, we report: the number of lines of code (in the C source), the total number of loops, the number of loops that are correctly bounded by each tool, and the computation time. We do not report the computation time for SWEET because we only had access to it through an online applet. For oRange, computation time is below the measurement resolution (10ms), except for edn, where it reaches 50ms. We ran PAGAI with the -d pk -t lw+pf options. For the PolyBench benchmarks, we did not succeed in running them with SWEET due to the online applet limitation. For the correlation benchmark, we did not succeed in running it with PAGAI, it terminates without giving any result.

The execution time of Polymalys is typically higher than that of PAGAI because we introduce more variables and constraints. We believe that we can reduce the gap with additional optimizations, however Polymalys will probably remain more costly , because it works at a lower level of abstraction.

Cases where tools fail to analyze some loop bounds are depicted in bold. There is only one benchmark for which Polymalys did not find a loop bound: for janne complex. The difficulty is that it contains complex loop index updates inside a if-then-else. On the contrary, there are several cases where Polymalys successfully estimates loops bounds, while the other tools fail to do so. Note that PAGAI does not specifically compute loop bounds, instead it computes loop invariants. We deduced loop bounds from these invariants.

Loop bounds examples

We further illustrate the differences between tool capabilities on some synthetic program examples. The off-by-one error (lines 5-6) may cause the array cell assignment (line 7) to overwrite the bound variable with the value 100. Polymalys correctly detects that the loop may iterate 100 times, while oRange and SWEET detect a maximum of 10 iterations. PAGAI also bounds to 10 iterations, but warns about a possible undefined behavior and unsafe result. Note that the bound depends on the stack variable allocation layout. In our experiments, the compiler allocates the bound variable next to the array. Such an information is much easier to analyze at the binary code level than at the source code level.

Example 13. The following example shows the benefits of a relational domain: for (i = base; i < end; i++); } Here, we do not know statically the value of end and base. However, due to the if statement (line 4), Polymalys introduces the constraint endbase ≤ 10. Thus, Polymalys bounds the loop correctly (10 iterations), while PAGAI, oRange and SWEET do not.

Example 14. Finally, we report analysis results for the motivational example of Figure 1. Polymalys correctly finds that the loop bound is equal to the maximum size of the UDP payload; PAGAI, oRange and SWEET fail to provide any bound.

Conclusion

In this paper we propose a novel technique for performing abstract interpretation of binary code using polyhedra. It consists in adding new variables to the polyhedra as the analysis progresses, and maintaining a correspondence with registers and memory addresses. Thanks to the relational properties of polyhedra, our technique naturally provides information on pointer relations when compared to other techniques based on non-relational domains. While the complexity of our method is currently still higher than other existing techniques, we believe that there is room for improvement. In particular, we are planning to extend our work with a modular procedure analysis and a data-structure analysis.

Fig. 1 .

 1 Fig. 1. Network-inspired benchmark

  l1 : I1, l2 : I2, . . . , ln : END Labels (L) ::= {l1, l2, . . .} Registers (R) ::= {r1, r2, . . .} Constants (C) ::= {c1, c2, . . .} Instructions (I)

Fig. 4 .

 4 Fig. 4. Semantics of control-flow operations.

r) end if 6: end for 7 :

 7 for all a ∈ Dom( * 1 ) do 8: if * 1 (a) = * 2 (a) then * (a) ← * 1 (a) end if 9: end for 10: return (p, R , * )

Example 11 . 1 #Author

 111 The following example contains pointer aliasing and pointer arithmetic: foo() { int i, bound = 10; int *ptr = &bound; ptr++; ptr--; *ptr = 15; k = 0; for (i = 0; i < bound; i++); } PAGAI does not find the loop bound (the loop is considered unbounded), because it does not infer that ptr = &bound when executing the instruction *ptr=15. Other tools bound the loop correctly (15 iterations). Example 12. The following example contains an off-by-one array access: Version of « Static Analysis Of Binary Code With Memory Indirections Using Polyhedra», Clément Ballabriga, Julien Forget, Laure Gonnord , Giuseppe Lipari, Jordy Ruiz accepted at VMCAI'19 int tab[SIZE]; if ((offset > SIZE) || (offset < 0)) return -1; tab[offset] = 100; for (i = 0; i < bound; i++); }

  end, i; if (end -base > MAXSIZE) end = base + MAXSIZE;

  The previous content at this address is replaced by the content of the source register (see Replace below). Otherwise, we create a new memory mapping (see Create below). An alias free state contains at most one address variable that must-alias with r 1 . It may however contain several mayalias address variables a . For each such a , this means that a either equals r 1 , which requires a Replace, or is different from r 1 , which has no impact. We apply operator on both cases to manage this uncertainty, and add the constraints for each may-alias address (see May below).

	(STORE r 1 r 2 ) =	λs.Replace(a)(May(s)) if ∃a ∈ vars A (p), a ≡ r 1 λs.Create(May(s)) otherwise
	With ( denotes function composition):

  (r) by R 1 (r) in p 2 , R 2 ,and * 2 13: end for 14: return (p 2 , R 2 , * 2 )

			R 2 , * 2 )
	2: Vc ← vars(p1) ∩ vars(p2)	common variables
	3: p ← proj (p1, Vc)	proj (p2, Vc)
	4: B ← npiv (p )
	5: for all (xi, xj) ∈ varsA(p1) × varsA(p2) do
	6:
			2
	9:	end if
	10: end for
	11: for all r ∈ Dom(R 1 ) ∩ Dom(R 2 ) do	variables are trivially equivalent
	12:	Replace R 2

vi = linexpr (xi, p1, B); vj = linexpr (xj, p2, B) 7: if vi = [] and vj = [] and vi = vj then variables are equivalent 8: Replace xj by xi and * (xj) by * (xi) in p 2 , R 2 , and *

Table 1 .

 1 Benchmark results.

				Loops Correctly Bounded	Time (ms)
	Benchmark	LoC Loops Polymalys SWEET PAGAI oRange Polymalys PAGAI
	crc	16	1	1	1	1	1	150	40
	fibcall	22	1	1	1	1	1	230	50
	janne complex 26	2	1	2	1	1	870	140
	expint	56	3	3	2	3	3	732	9140
	matmult	84	5	5	5	5	5	3455	1380
	fdct	149	2	2	2	2	2	7421	2150
	jfdctint	165	3	3	3	3	3	10660	1960
	fir	189	2	2	2	2	1	4989	390
	edn	198	12	12	12	9	12	21356	15660
	ns	414	4	4	4	4	4	1700	380
	gemver	186	10	10	N/A	10	10	12136	6029
	covariance	138	11	11	N/A	11	11	7248	836
	correlation	168	13	13	N/A N/A	13	9129	25062
	nussinov	143	8	8	N/A	8	8	7272	2811
	floyd-warshall 112	7	7	N/A	2	7	2904	468

Partially funded by the French National Research Agency, Corteva project (ANR-17-CE25-0003), and CODAS project (ANR-17-CE23-0004-01)

Author Version of « Static Analysis Of Binary Code With Memory Indirections Using Polyhedra», Clément Ballabriga, Julien Forget, Laure Gonnord , Giuseppe Lipari, Jordy Ruiz accepted at VMCAI' [START_REF] Hind | Pointer analysis: Haven't we solved this problem yet?[END_REF] Author Version of « Static Analysis Of Binary Code With Memory Indirections Using Polyhedra», Clément Ballabriga, Julien Forget, Laure Gonnord , Giuseppe Lipari, Jordy Ruiz accepted at VMCAI'19 Lemma 5. Let a 1 , a 2 ∈ A. We have: γ(a 1 ) ∪ γ(a 2 ) ⊆ γ(a 1 a 2 ).

Proof. From Lemma 1 and Lemma 4. Lemma 6. Let a 1 , a 2 ∈ A. We have: γ(a 1 ) ∪ γ(a 2 ) ⊆ γ(a 1 a 2 ) Proof. From Lemma 5, Lemma 4 and Property 3.

Lemma 7. Let a ∈ A. We have: γ(a) ⊆ γ(antialias(a))

The soundness of antialias follows.

Lemma 8. Let P be a MEMP program. Let M = Interpret(P ). Then, for all labels , of P :

Proof. Trivially follows from the formal semantics and from the definition of transfer functions, except for STORE. Let a = (p , R , * ) = (STORE r 1 r 2 ) (a). The proof follows from noting that: 1) Both in the Create and Replace cases, we obtain * (R (r 1 )) = R (r 2 ), which is coherent with the formal semantics of STORE;

2) The soundness of May follows from the soundness of and Replace.

Lemma 9. Algorithm 1 terminates.

Proof. Because is applied on loop headers and is a valid widening operator.

Related works

Abstract interpretation using polyhedra has been first described in [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF]. Static analysis tools such as Astree [START_REF] Kästner | Astrée: Proving the absence of runtime errors[END_REF], Frama-C [START_REF] Correnson | Combining analyses for c program verification[END_REF] or PAGAI [START_REF] Henry | Pagai: A path sensitive static analyser[END_REF] use various abstract domains (including polyhedra) to generate invariants for proving various properties, such as the absence of array out-of-bounds accesses for instance. While Astree and Frama-C work on the Abstract Syntax Tree, PAGAI processes LLVM Intermediate Representation (IR). Compared to our approach, both the AST and LLVM representations are closer to the source code, and contain information on variables and their types, and also a precise control flow. This makes the analysis easier to design, but less precise as far as WCET is concerned.