A Ragouilliaux

G Ovarlez

N Shahidzadeh-Bonn

Benjamin Herzhaft

T Palermo

P Coussot

Transition from a simple yield stress fluid to a thixotropic material

Keywords: 82.70.-y, 83.60.Pq, 61.43.Hv

From MRI rheometry we show that a pure emulsion can be turned from a simple yield stress fluid to a thixotropic material by adding a small fraction of colloidal particles. The two fluids have the same behavior in the liquid regime but the loaded emulsion exhibits a critical shear rate below which no steady flows can be observed. For a stress below the yield stress, the pure emulsion abruptly stops flowing, whereas the viscosity of the loaded emulsion continuously increases in time, which leads to an apparent flow stoppage. This phenomenon can be very well represented by a model assuming a progressive increase of the number of droplet links via colloidal particles.

I. Introduction

Various systems such as foams, emulsions or colloidal suspensions, are made of elements interacting via soft interactions, i.e. which vary progressively with the distance between the elements [START_REF] Coussot | Rheometry of pastes, suspensions and granular materials[END_REF][START_REF] Trappe | [END_REF]. Their jammed character, namely their ability to undergo a solid-liquid transition when submitted to a sufficient stress, arouses the interest of physicists as it appears typical of a fourth state of matter with some analogy with glass behavior [3]. Their rheological behavior is usually described with the help of simple yield stress models, which may yield velocity profiles with unsheared regions when the stress distribution is heterogeneous. Recently some works showed that the flow of such materials exhibits "true shear-banding", i.e. two regions with significantly different shear rates in geometries in which the shear stress is homogeneous [4]. Such a behavior is not predicted by usual yield stress models for which the transition between the solid and the liquid phases is smooth.

Two questions thus emerge: (i) how can a material become a "shear-banding" material? and (ii) how does the transition towards shear-banding occur in time for a given material? The first question was recently addressed by Bécu et al. [5] who showed that a simple yielding emulsion can exhibit shear-banding if it is made "adhesive", i.e. the droplets develop some adhesive interactions. Here we focus on question (ii), namely the time-dependency of this shear-banding effect and its origin, as it is likely that the transition from the purely solid (starting from rest) or purely liquid (starting from an intense preshear) regimes towards shearbanding is a complex phenomenon from which useful information for our understanding of these materials could be extracted. Remark that there is certainly some interplay of these phenomena with the thixotropic character of these fluids (time variations of the apparent viscosity). For example it was shown that a clay-water system can exhibit very different behavior types, and in particular shear-banding or even fractures, depending on the time it spent at rest before undergoing a stress [6].

Here we compare the behavior of an emulsion either pure or loaded with colloidal particles. In order to have a straightforward appreciation of the effective behavior of these materials we observe their flow characteristics in time with the help of MR (Magnetic Resonance) velocimetry. We show that the pure emulsion can be considered as a simple yield stress fluid with negligible thixotropy while the loaded emulsion is significantly thixotropic. Both fluids have the same behavior in the liquid regime and the thixotropy of the loaded emulsion is due to a progressive aggregation of droplets via colloidal particle links, which tends to increase the viscosity of the material which finally stops. As a result the apparent yield stress (below which no steady flow can be obtained) of the loaded emulsion is larger than the pure emulsion. We finally show that during stoppage the shear rate exhibits a specific evolution which is consistent with the predictions of a model assuming a progressive aggregation of droplets in time.

II. Materials and procedures

We prepared the pure emulsion by progressively adding the water in an oil-surfactant solution (Sorbitan monooleate, 2%) under high shear. The size distribution of water droplets was almost uniform around a diameter of 1 μm % 20 ±

. The water droplet concentration was 70%, a concentration larger than the maximum packing fraction of uniform spheres, which means that the droplets were probably slightly deformed. However the droplet deformation due to flow can be neglected during all the flows considered here since, due to the very small droplet size, the interfacial stress is much larger than the stress due to flow. As a consequence the average shape of droplets can be considered as constant during all our tests. outer cylinder radius: 6cm; height of sheared fluid: 11cm) inserted in a MRI set-up. We imposed the rotation velocity of the inner cylinder and measured the torque and the velocity profile ( ) (r V in which r is the distance from the central axis) inside the sample at successive times by magnetic resonance velocimetry. For a given rotation velocity the torque ( M ) was constant in time even if for the loaded emulsion the velocity profile significantly evolved (see below). The principle of the magnetic resonance velocimetry technique used was described elsewhere [7]. The typical time needed to get one velocity profile was 15s. The uncertainty on the local velocity measurements was typically We also carried out NMR density tests with the help of a new technique [8] for separating the oil and the water signals, which make it possible to detect within 1% some heterogeneity of the water fraction resulting from migration effects. No migration was observed in the pure emulsion. With the clayey emulsion no such migration effect was observed during the first 20min. of flow but a significant migration developed beyond 30min. of flow. In the following we present only the data corresponding to the first stage (with negligible migration). We also ) with rough surfaces (roughness: 0.18mm).

carried

III. Experimental results

A. Solid regime

From creep tests we find that both materials exhibit a solid and a liquid regimes: at low stress the induced deformation tends to saturates (solid regime), whereas at high stress they tend to flow at a constant rate (liquid regime). This is the typical behavior observed for pasty materials [9]. Let us first compare the material properties in their solid regime from dynamic tests (frequency: 1Hz, strain amplitude: 1%). For the pure emulsion the elastic modulus ( ' G )

slightly increases with time (see inset of Fig. 1) (i.e. of 7% after 3000s), which means that, as a first approximate, it negligibly ages in its solid regime. Thus the possible evolution of the droplet configuration in time only induces a slight increase of the network strength. In contrast for the loaded emulsion ' G significantly increases in time (i.e. of 70% after 3000s).

Thus, as time goes on, more and more clay links form between droplets, which strengthens the network of droplet contacts. ) for the pure emulsion. Inset: Elastic modulus (scaled by the value reached at 40s) of the pure (circles) and clayey (squares) emulsions as a function of time.

B. Liquid regime

Let us now consider the behavior in the liquid regime. In that aim we preshear the material at a large velocity and leave it at rest few seconds, then we impose a given rotation velocity and observe the velocity profile inside the gap. During a short initial period (say less than few seconds) the rotation velocity varies as a result of inertia and solid-liquid transition effects.

Afterwards, for the pure emulsion the velocity profile is constant. This means that this material exhibits negligible thixotropy in the liquid regime: the specific (average) droplet configuration associated with a given flow rate is reached almost immediately. In contrast, with the loaded emulsion the velocity profile significantly evolves in time: it progressively moves towards the central axis (cf. Fig. 2) and finally reaches a steady shape. During this process the outer region apparently progressively stops flowing but there still might remain some slow flow at velocities smaller than the uncertainty on these measurements.

In a Couette geometry the shear stress distribution (τ as a function of r ) in the material is well known from the momentum equation and we have ) are not represented, which explains why there is no data beyond some critical distance for each profile.

A. Data analysis

Let us now compare the constitutive equation of these materials in the liquid regime in steady state as it can be deduced from a set of experiments under different torques from MRI data, by computing the local shear rate and shear stress for each value of r as above described. The flow curve (steady state τ vs γ ) of the pure emulsion is that of a simple yield stress fluid, i.e. the shear stress tends to a plateau at low shear rates below which no steady flows are obtained and a Herschel-Bulkley model very well represents the flow curve (see Fig. 1). The loaded emulsion exhibits a behavior similar to that of the pure emulsion for shear rates larger than say -1 s 10 . This means that the clay links no longer play any role and suggest that there could be all broken.

The striking difference between the two materials concerns the minimum stress level allowing steady state flow: the yield stress of the clayey emulsion is significantly larger than that for the pure emulsion and just beyond this yield stress the fluid flows at a shear rate larger than c γ . Forgetting the short plateau of the flow curve for the clayey emulsion we get the following rough picture: both materials have a similar behavior in the liquid regime but the flow curve of the clayey emulsion is truncated at a critical shear stress larger than the yield stress of the pure emulsion, and no stable flows can be obtained at a shear rate below the critical value associated with this yield stress. For such materials shear-banding will occur in shear flows at apparent shear rates smaller than c γ . This behavior is reminiscent of the trend observed by Becu et al. [5] with an adhesive emulsion, but here we get a precise picture of the different fluid behavior types at a local scale and we show that the difference finds its origin in some time-dependent effect.

With regards to our scheme of the microstructure it is natural to consider that this behavior is the result of a dynamic effect: a progressive transition from flow to stoppage, which is likely due to the progressive formation of larger aggregates. Thus as a first approximation we can consider that all the clay links are broken during a rapid shear, so that the clay fraction plays a negligible role (in agreement with our observation of the similar behavior of the two materials in the liquid regime), but below a critical stress the number of links increases significantly, which slows down the flow so that the fraction of links further increases and so on until full stoppage.

Previous studies of the thixotropic properties in the liquid regime generally relied on macroscopic data under poorly controlled conditions. Here, from MR velocimetry we have straightforward information concerning the local viscosity variations in time (see Fig. 3). We can remark that the shear rate vs time curves for different stress values below c τ look similar in a logarithmic scale: they start with a plateau at short time then rapidly decrease after some time, which decreases as the stress decreases. This similarity definitely appears when plotting the shear rate vs the time scaled by appropriate factors, respectively ) ( 0 τ γ and ) (τ θ (see Fig. 4). This is reminiscent of an effect recently observed [11]: for various suspensions the elastic modulus vs time curves after the liquid-solid transition under different stress fall along a master curve when the elastic modulus and the time are scaled by appropriate functions of the stress. Here we get the following general result: for a stress imposed below the critical stress the shear rate evolves as

) ( 0 θ γ γ t f = , with [ ] c c γ τ τ τ τ γ ) ( ) ( 0 0 0 - - ≈ and
) ( 1 τ τ θ - ≈ c k
. This similarity of the behavior strongly suggests that the structural state of the material is a function of

θ ξ t =
only, which thus can be considered as describing the physical age, i.e. the state of structure, of the material during such a test. It is remarkable that the evolutions of this parameter only depend on the stress. This contrasts with usual models (e.g. [6]) in that field which assume that the structure parameter essentially depends on the shear rate.

From the above equations we deduce the apparent behavior of the material for all ξ :

( ) ) ( ) ( 0 0 ξ γ ξ η τ τ f + = (1)
where

c c γ τ τ η ) ( 0 0 - = .
In this equation the first term, i.e. 0 τ , is independent of the structural state and thus can be considered as reflecting the yield stress of the clayey emulsion when there is no droplet-clay links. This yield stress takes its origin in the jammed droplets network with clay possibly affecting the interfacial tension between oil and water, which explains that it may differ from the yield stress of the pure emulsion. The second term (the "viscous term") is associated with both the viscous effects and the droplet-clay links in the emulsion. It is worth noting that the behavior described by equation [START_REF] Coussot | Rheometry of pastes, suspensions and granular materials[END_REF] does not correspond to a Bingham model. Indeed the factor of the shear rate ( ) ( 0 ξ η f ), i.e. the so-called "plastic viscosity" which is constant in the Bingham model, in our model depends on the structural state ξ and thus depends on the shear stress.

B. Model

The above analysis suggests that one may see the clayey emulsion as a suspension of linked elements in an interstitial homogeneous yield stress fluid. We will make two simplifying assumptions: (i) the stress in this suspension may be expressed as the sum of a "yielding" term strictly due to the jamming-unjamming process during flow and a simple viscous term with an apparent viscosity (η ), resulting from the additional viscous dissipation due to flow and function of the actual concentration of aggregates; (ii) the presence of aggregates negligibly affects the yielding part. Such a picture is qualitatively in perfect agreement with the conclusions of our data interpretation (see above). In this context the viscous term may be found from the viscosity of a Newtonian suspension in a liquid of viscosity 0 η .

We assume that at a given time the droplet-clay links form more or less large (rigid) networks which trap some volume of the interstitial liquid. Various models exist for the viscosity of rigid spherical inclusions [12] but it seems more relevant to use the expression found from energy dissipation considerations in the context of purely viscous suspensions of aggregative particles [13]: 4). Remark that the aggregate size tends towards infinity at stoppage, so that the above theory no longer holds and the corresponding stress contributes to increase the yield stress of the material, which explains that the apparent yield stress of the clayey emulsion is larger than pure emulsion. The above equation provides the relation between the effective advancement of the aggregation (via agg φ ) and the time and the stress applied. It is remarkable that the stress is solely involved in this relationship whereas in general for such physical phenomena the characteristic time is proportional to γ 1 . One finds a possible beginning of explanation of such an effect in the usual description of simple liquids [START_REF] Glasstone | The Theory of Rate Processes[END_REF], glassy [START_REF] Monthus | [END_REF] or jammed systems [16] in which the elements are assumed to be in a potential well with an energy barrier of depth E Δ . We can assume that such a frame applies for the clay particles involved in claydroplet links. The probability of escape from this potential well due to thermal agitation is

( ) 2 0 ) ( 1 ) 1 (
T k E B Δ - exp
per attempt, in which B k is the Boltzmann constant and T the temperature. In this context the influence of stress is to decrease the energy barrier by a factor d τ in which d is a characteristic distance of the system, and thus to foster the fluid flow. In order to have the possibility to form a clay bridge two droplets must be neighbours. Over a time of flow t Δ the number of occurrence of droplet neighbouring is proportional to γ but the time during which two droplets can be considered as neighbours is inversely proportional to γ , so that the rate of flow has no influence on the number of links (this result is valid only because we neglected possible link breakages). To sum up we find that the rate of variation of the number of links increases with the stress and is independent of the shear rate, in agreement with our experimental results. However a quantitative prediction of the effect of stress on θ requires a more complete theoretical approach.

V. Conclusion

We have shown that the transition from a non-thixotropic to a thixotropic emulsion occurs by the disappearance of steady flows below a critical stress associated with a critical shear rate.

Here this process is well explained by the increase of droplet links via clay particles. These effects might be general for jammed systems for which the aging is due to reversible particle aggregation.

5 -

 5 m/s. Slight wall slip effects possibly occurred but did not affect the transient measurements as they remained constant for a given torque and did not affect our rheological interpretations based on measurements inside the sample.

Figure 1 :

 1 Figure 1: Steady-state shear stress vs shear rate for the pure (circles) and clayey (squares) emulsions as deduced from local measurements within the wide-gap rheometer under different torque values. The solid line is a Herschel-Bulkley model ( n e Kγ

.

  e. at a given torque level the stress decreases with the distance from the axis. Our observations concerning the velocity profiles under a constant torque imply that in the regions of low stresses (near the outer cylinder) the material tends to stop flowing while in the regions of large stresses (near the inner cylinder) it tends to reach a steady-state flow. This effect is reminiscent of the viscosity bifurcation effect observed with various material types[10]. However here, in contrast with previous studies, we have a straightforward information concerning the effective (local) shear rate in time. Indeed we can compute the local shear rate (γ ) for any stress (associated with a given distance r ) at any time: we have Looking now at the evolutions in time of the shear rate under different stress values (i.e. at different distances) we effectively observe the viscosity bifurcation effect (see Fig.3): for c τ τ > the shear rate tends to a constant value, for c τ τ < the shear rate progressively tends to zero, with a characteristic time increasing with τ . As a consequence there exists a critical shear rate ( c γ ), associated with c τ , below which no steady-state flows can be obtained, the material tends to stop. From our data we get

Figure 2 :

 2 Figure 2: Velocity profiles at different times in the wide-gap Couette geometry as measured by MR velocimetry for the clayey emulsion for a rotation velocity of 15rpm: (from right to left) 30, 60, 90, 120, 150, 180, 210, 240, 270, 330, 390, 450, 540, 660, 780, 900s. The velocity values below the uncertainty ( -1 5 m.s 10 . 5 -) are not represented, which explains why there is no data beyond some critical distance for each profile.

Figure 3 :

 3 Figure 3: Effective shear rate as a function of time for different local stresses as determined from the MR velocity profiles for the clayey emulsion for a rotation velocity of 15rpm: (from top to bottom) 34.2, 31.7, 30.1, 29.4, 28.7, 28, 27.3, 26.7, 26.1, 25.5, 24.9, 24.3, 23.8, 22.8, 22.3, 21.8, 21.4, 20.9, 20.5, 20.1, 19.7, 18.9 Pa.

  φ is a parameter which represents the volume fraction of rigid inclusions formed by interstitial liquid trapped in clay-droplet networks. The specific maximum packing fraction m φ used in the Mills model[13] is equal to 0.57. A different choice for this value would slightly change the value of the fitting parameter in the following analysis but not the qualitative results. This viscous term can be associated with both terms describe the apparent viscosity of the suspension in a simple liquid of viscosity 0 η , which provides a correspondence between agg φ and ξ . For example a kinetics relation of the form to get a good fit of the above equation to the master curve (see Fig.

Figure 4 :

 4 Figure 4: Shear rate (data of Figure 3 for Pa 3 . 27 ≤ τ


Claude-Pierre Jeannerod
email: claude-pierre.jeannerod@ens-lyon.fr

Herv

Christophe Knochel
email: christophe.monat@st.com

Guillaume Monat

Revy
email: grevy@eecs.berkeley.edu

C.-P Jeannerod

• H Knochel

C Monat

G Revy

Computing floating-point square roots via bivariate polynomial evaluation

Keywords: Binary floating-point arithmetic, square root, correct rounding, IEEE 754, polynomial evaluation, instruction-level parallelism, rounding error analysis, C software implementation, VLIW integer processor

In this paper we show how to reduce the computation of correctly-rounded square roots of binary floating-point data to the fixed-point evaluation of some particular integer polynomials in two variables. By designing parallel and accurate evaluation schemes for such bivariate polynomials, we show further that this approach allows for high instruction-level parallelism (ILP) exposure, and thus potentially low latency implementations. Then, as an illustration, we detail a C implementation of our method in the case of IEEE 754-2008 binary32 floating-point data (formerly called single precision in the 1985 version of the IEEE 754 standard). This software implementation, which assumes 32-bit unsigned integer arithmetic only, is almost complete in the sense that it supports special operands, subnormal numbers, and all rounding-direction attributes, but not exception handling (that is, status flags are not set). Finally we have carried out experiments with this implementation on the ST231, an integer processor from the STMicroelectronics' ST200 family, using the ST200 family VLIW compiler. The results obtained demonstrate the practical interest of our approach in that context: for all rounding-direction attributes, the generated assembly code is optimally scheduled and has indeed low latency (23 cycles).

INTRODUCTION

T HIS paper deals with the design and software im- plementation of an efficient sqrt operator for computing square roots of binary32 floating-point data. As mandated by the IEEE 754 standard (whose initial 1985 version [START_REF]IEEE standard for binary floating-point arithmetic[END_REF] has been revised in 2008 [START_REF]IEEE standard for floating-point arithmetic[END_REF]), our implementation supports gradual underflow and the four roundingdirection attributes required for binary format implementations. However, the status flags used for handling exceptions are not set.

As for other basic arithmetic operators, the IEEE 754 standard specifies that sqrt operates on and returns floating-point data. Floating-point data are either special data (signed infinities, signed zeros, not-a-numbers (NaN)) or nonzero finite floating-point numbers. In radix two, nonzero finite floating-point numbers have the form x = ± m • 2 e , with e an integer such that

e min ≤ e ≤ e max , (1)
and m a positive rational number having binary expansion m = (0.0

• • • 0 λ zeros 1m λ+1 • • • m p-1 ) 2 . ( 2
)
Since m is nonzero, the number λ of its leading zeros varies in {0, 1, . . . , p-1}. If |x| < 2 emin then x is subnormal else it is normal. On the one hand, subnormal numbers are such that e = e min and λ > 0. On the other hand, normal numbers will be considered only through their normalized representation, that is, the unique representation of the form ±m • 2 e for which λ = 0.

The parameters e min , e max , p used so far represent the extremal exponents and the precision of a given binary format. In this paper, we assume they satisfy e min = 1 -e max and 2 ≤ p ≤ e max .

This assumption is verified for all the binary formats defined in the IEEE 754-2008 standard [START_REF]IEEE standard for floating-point arithmetic[END_REF]. This standard further prescribes that the operator sqrt : x → r specifically behaves as follows:

• If x equals either -0, +0, or +∞ then r equals x.

• If x is nonzero negative or NaN then r is NaN. Those two cases cover what we shall call special operands. In all other cases x is a positive nonzero finite floatingpoint number, that is,

x = m • 2 e , (3)
with e as in [START_REF]IEEE standard for binary floating-point arithmetic[END_REF] and m as in [START_REF]IEEE standard for floating-point arithmetic[END_REF]; the result specified by the IEEE 754-2008 standard is then the so-called correctlyrounded value

r = •( √ x), (4)
where • is any of the four rounding-direction attributes: to nearest even (RN), up (RU), down (RD), and to zero (RZ). In fact, since √ x ≥ 0, rounding to zero is the same as rounding down. Therefore considering only the first three rounding-direction attributes is enough:

• ∈ {RN, RD, RU}.
As we will see in Section 2, deducing r from x essentially amounts to taking the square root, up to scaling, of the significand m. For doing this, many algorithms are available (see for example the survey [START_REF] Montuschi | Survey of square rooting algorithms[END_REF] and the reference books [START_REF] Markstein | IA-64 and Elementary Functions: Speed and Precision[END_REF], [START_REF] Cornea | Scientific Computing on Itanium R -based Systems[END_REF], [START_REF] Ercegovac | Digital Arithmetic[END_REF]).

The method we introduce in this paper is based exclusively on fixed-point evaluation of a suitable bivariate polynomial with integer coefficients. Since polynomial evaluation is intrinsically parallel, this approach allows for very high ILP exposure. Thus, in some contexts such as VLIW integer processors, a significant reduction of latency can be expected. For the ST231, a STMicroelectronics VLIW processor which does not have native floating-point capabilities, we show in this paper that this is indeed the case: the assembly codes generated for optimized implementations of our approach turn out to be optimally scheduled and have latencies reduced by over 30% compared to previously fastest available methods, implemented in [START_REF] Jeannerod | Faster floating-point square root for integer processors[END_REF]; also, the latency overhead for subnormal support is only 2 cycles, yielding a latency of 23 cycles for all •. This latency for square root compares favourably with the latencies of addition and multiplication which, in the same context, range respectively from 23 to 27 and 18 to 21 cycles, depending on • (see http://flip.gforge.inria.fr/ and [8, p. 4]).

The paper is organized as follows. Section 2.1 details three mathematical properties of square roots of binary floating-point numbers, and Section 2.2 shows how to use them to deduce the usual high level algorithmic description of the sqrt operator.

In Section 3.1 we then show how to introduce suitable bivariate polynomials that approximate our square root function. In particular, we give some approximation and evaluation error bounds that are sufficient to ensure correct rounding, along with an example of such a polynomial and its error bounds in the case of the binary32 format. Section 3.2 then details, for each roundingdirection attribute, how to deduce a correctly-rounded value from the approximate polynomial value obtained so far. A summary of our new approach is given in Section 3.3.

A standard C99 implementation of this approach is given in Section 4, for the binary32 format and assuming that 32-bit unsigned integer arithmetic is available: Section 4.1 shows how to handle special operands; Section 4.2 deals with the computation of the result exponent and a parity bit related to the input exponent (which is needed several times in the rest of the algorithm); Sections 4.3 and 4.4 show how to compute the evaluation point and the value of the polynomial at this evaluation point; there, we also explain how the accuracy of the evaluation scheme has been verified; Finally, Section 4.5 details how to implement correct rounding. This results in three separate square root implementations, one for each rounding-direction attribute.

For these implementations, all that is needed is a C compiler that implements 32-bit arithmetic. However, our design has been guided by a specific target, the ST231 VLIW integer processor from STMicroelectronics' ST200 family. Section 5 is devoted to some experiments carried out with this target and the ST200 family VLIW compiler. After a review of the main features of the ST231 in Section 5.1, the performance results we have obtained in this context are presented and analysed in Section 5.2.

PROPERTIES OF FLOATING-POINT SQUARE

ROOTS AND GENERAL ALGORITHM

The general scheme for moving from (3) to (4) essentially follows from three properties of the square root function in binary floating-point arithmetic.

Floating-point square root properties

Property 2.1: For x in (3), the real number √ x lies in the range of positive normal floating-point numbers, that is,

√ x ∈ 2 emin , (2 -2 1-p ) • 2 emax .
This first property (see [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF] for a proof) implies that •( √ x) is a positive normal floating-point number. Correctly-rounded square roots thus never denormalize nor under/overflow, a fact which will simplify the implementation considerably.

In order to find the normalized representation of

•( √ x), let e = e -λ and m = m • 2 λ . (5)
It follows that the positive (sub)normal number x defined in (3) satisfies x = m • 2 e and that m ∈ [1, 2).

Taking the square root then yields

√ x = • 2 d ,
where the real and the integer d are given by

= σ √ m with σ = 1, if e is even, √ 2, if e is odd, (6)
and, using to denote the usual floor function,

d = e /2 . (7)
It follows from the definition of σ and m that is a real number in [START_REF]IEEE standard for binary floating-point arithmetic[END_REF][START_REF]IEEE standard for floating-point arithmetic[END_REF]. Therefore, rounding √ x amounts to round , and we have shown the following property: Property 2.2: Let x, , d be as above. Then

•( √ x) = •( ) • 2 d .
In general the fact that ∈ [1, 2) only implies the weaker enclosure •( ) ∈ [START_REF]IEEE standard for binary floating-point arithmetic[END_REF][START_REF]IEEE standard for floating-point arithmetic[END_REF]. This yields two separate cases: if •( ) < 2 then Property 2.2 already gives the normalized representation of the result r; if •( ) = 2 then we must further correct d after rounding, in order to return r = 1 • 2 d+1 instead of the unnormalized representation 2 • 2 d given by Property 2.2. The next property, proven in [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF], characterizes precisely when such a correction is needed or not. In particular, it is never needed in rounding-to-nearest mode.

Property 2.3: One has •( ) = 2 if and only if • = RU and e is odd and m = 2 -2 1-p .

High level description of square root algorithms

Together with the IEEE 754-2008 specification recalled in Introduction, the properties of Section 2.1 lead to a general algorithm for computing binary floating-point square roots, shown in Figure 1 and Figure 2 for the different rounding-direction attributes. In particular, •( ) and d are two functions of m and e which can be computed independently from each other. (We will see in Section 4.2 that the correction step in Figure 2 can in fact be avoided when using the standard encoding of binary floating-point data.)

x = m • 2 e Handle special input Compute d Compute •( ) x ∈ {±0, ±∞, NaN, x < 0} or r = •( ) • 2 d r ∈ {±0, +∞, NaN} or
Given m and e, computing d is algorithmically easy since by ( 5) and [START_REF] Jeannerod | Faster floating-point square root for integer processors[END_REF] we have

d = (e -λ)/2 . (8)
However, computing •( ) from m and e is far less immediate and typically dominates the cost. In the next section, we present a new way of producing •( ), which we have chosen because we believe it allows to express the most ILP.

COMPUTING •( ) BY CORRECTING TRUN-CATED APPROXIMATIONS

It is useful to start by characterizing the meaning of •( ) for each rounding-direction attribute •. Since the real number belongs to [1, 2), we have the following, where n is a normal number:

• RN( ) is the unique n such that -2 -p < -n < 2 -p , (9)
• RD( ) is the unique n such that

0 ≤ -n < 2 1-p , (10)
• RU( ) is the unique n such that

-2 1-p < -n ≤ 0. (11)
Both inequalities in (9) are strict, since the square root of a floating-point number cannot be exactly halfway between two consecutive floating-point numbers (see [START_REF] Markstein | IA-64 and Elementary Functions: Speed and Precision[END_REF]Theorem 9.4], [5, p. 242], or [6, p. 463]). Note also that ( 9) and ( 10) both imply n ∈ [1, 2), while [START_REF] Pi | High-speed double-precision computation of reciprocal, division, square root and inverse square root[END_REF] implies the weaker enclosure n ∈ [START_REF]IEEE standard for binary floating-point arithmetic[END_REF][START_REF]IEEE standard for floating-point arithmetic[END_REF]. Given p, •, m , and σ, there are many ways of producing the bits of such an n. A way that will allow to express much ILP is by correcting a truncated approximation of . This approach (detailed for example in [6, p. 459] for division) has three main steps:

• compute a "real number" v approximating from above with absolute error less than 2 -p , that is,

-2 -p < -v ≤ 0; (12)
• deduce u by truncating v after p fraction bits:

0 ≤ v -u < 2 -p ; (13)
• obtain n by adding, if necessary, a small correction to u and then by truncating after p -1 fraction bits. Of course the binary expansion of v in (12) will be finite: by "real number" we simply mean a number with precision higher than the target precision p. Typically, v will be representable with at most k bits, with k the register width (for example, p = 24 and k = 32 for our implementation in Section 4). On the other hand, using

≤ √ 2 • √ 2 -2 1-p one may check that if v satisfies (12) then necessarily 1 ≤ v < 2. ( 14
)
Therefore, the binary expansion of v has the form

v = (1.v 1 . . . v p-1 v p . . . v k-1 ) 2 . (15)
Our approach for computing v as in ( 12) and ( 15) by means of bivariate polynomial evaluation is detailed in Section 3.1 below. Once v is known, truncation after p fraction bits gives

u = (1.v 1 . . . v p-1 v p ) 2 . ( 16
)
The fraction of u is wider than that of n by one bit. We will see in Section 3.2 how to correct u into n by using both this extra bit and the fact that, because of ( 12) and ( 13),

| -u| < 2 -p . ( 17
)

Computing v by bivariate polynomial evaluation

The problem is, given p, m , and σ, to compute an approximation v to such that (12) holds. Such a value v will in fact be obtained as a solution to

( + 2 -p-1 ) -v < 2 -p-1 . (18)
Although slightly more strict than [START_REF] Raina | FLIP: a floating-point library for integer processors[END_REF], this form will be the natural result of some derivations based on the triangle inequality.

Computing v such that [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF] holds usually relies on iterative refinement (Newton-Raphson or Goldschmidt method [6, §7.3]), or univariate polynomial evaluation [START_REF] Jeannerod | Faster floating-point square root for integer processors[END_REF] (see also [START_REF] Agarwal | Series approximation methods for divide and square root in the Power3 TM processor[END_REF] for a different context, namely when an FMA instruction is available), or a combination of both [START_REF] Pi | High-speed double-precision computation of reciprocal, division, square root and inverse square root[END_REF], [START_REF] Raina | FLIP: a floating-point library for integer processors[END_REF]. The approach we shall detail now is based exclusively on polynomial evaluation. However, instead of using two polynomials as in [START_REF] Jeannerod | Faster floating-point square root for integer processors[END_REF], we will use a single one, which is bivariate; this makes the approach simpler, more flexible, and eventually faster provided some parallelism is available.

The main steps for producing v via bivariate polynomial evaluation are shown on the diagram below:

+ 2 -p-1 ----→ F (s, t) ?     function approximation v ← ------------- polynomial evaluation P (s, t)
First, using ( 6) and defining

τ = m -1, (19)
the real number + 2 -p-1 can be seen as the value at (s, t) = (σ, τ ) of the bivariate function

F (s, t) = 2 -p-1 + s √ 1 + t. (20)
Then, let

S = {1, √ 2} and T = {h • 2 1-p } h=0,1,...,2 p-1 -1
be the variation domains of σ and τ , respectively, and let

1 -= 1 -2 1-p .
Since T ⊂ [0, 1 -], a second step is the approximation of F (s, t) on {1, √ 2} × [0, 1 -] by a bivariate polynomial P (s, t). The function F being linear with respect to its first variable, a natural choice for P is

P (s, t) = 2 -p-1 + s • a(t), (21)
with a(t) a univariate polynomial that approximates

√ 1 + t on [0, 1 -].
The third and last step is the evaluation of P at (σ, τ ) by a finite-precision, straight-line program P, the resulting value P(σ, τ ) being assigned to v.

Intuitively, if a(t) is "close enough" to √ 1 + t over the whole interval [0, 1 -] and if P evaluates P "accurately enough" on the whole domain S × T then the resulting value v should be close enough to + 2 -p-1 in the sense of [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF]. This claim is made precise by the next lemma.

Lemma 1: Given p, σ, τ , a, P , P as above, let α(a) be the approximation error defined by

α(a) = max t∈[0,1 -] √ 1 + t -a(t) ,
and let ρ(P) be the rounding error defined by

ρ(P) = max (s,t)∈S×T |P (s, t) -P(s, t)| . Let further v = P(σ, τ ). If √ 2 • α(a) + ρ(P) < 2 -p-1 (22)
then v satisfies [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF]. Proof: Using the definitions of F and P , we have

( + 2 -p-1 ) -v = |F (σ, τ ) -P(σ, τ )| ≤ |F (σ, τ ) -P (σ, τ )| + |P (σ, τ ) -P(σ, τ )| ≤ σ √ 1 + τ -a(τ ) + ρ(P). Since 1 ≤ σ ≤ √ 2 and τ ∈ [0, 1 -], it follows that ( + 2 -p-1 ) -v ≤ √ 2 • α(a) + ρ(P),
which concludes the proof.

It remains to find a polynomial approximant a together with an evaluation program P so that α(a) and ρ(P) satisfy [START_REF] De Dinechin | Assisted verification of elementary functions using Gappa[END_REF]. In practice, since α(a) and ρ(P) may be real numbers, a certificate will consist in computing two dyadic numbers d α and d ρ such that

α(a) ≤ d α , ρ(P) ≤ d ρ , √ 2 • d α + d ρ < 2 -p-1 .
The construction of a and P is highly context-dependent: it is guided by both the value of p and some features of the target processor (register precision k, instruction set, latencies, degree of parallelism,...). The two paragraphs below illustrate how to choose a and P in the case where k = 32 and p = 24.

Constructing a polynomial approximant

Since P in [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF] will be evaluated at run-time, a small degree for a is usually preferred. One may guess the smallest possible value of deg(a) as follows. The rounding error ρ(P) in ( 22) being non-negative, a must satisfy

α(a) < 2 -p-3/2 . ( 23
)
Now let R[t] d be the set of univariate real polynomials of degree at most d and recall, for example from [13, §3.2], that the minimax degree-d approximation of

√ 1 + t on [0, 1 -] is the unique a * ∈ R[t] d such that α * d := α(a * ) ≤ α(a), for all a ∈ R[t] d . (24)
Thus, by combining ( 23) and ( 24), one must have

α * deg(a) < 2 -p-3/2 .
A lower bound on deg(a) can then be guessed by estimating α * d numerically for increasing values of d until 2 -p-3/2 is reached.

For example, for p = 24, the degree of a must satisfy α * deg(a) < 2 -25.5 . Comparing to the estimations 1 of α * d displayed in Table 1 indicates that a should be of degree at least 8.

Once we have an idea of the smallest possible degree for a, it remains to find a machine representation of a. For this representation, a typical choice is the monomial basis 1, t, t 2 , . . . together with some coefficients a 0 , a 1 , a 2 , . . . that are exactly representable using at most k bits.

1. Such estimations can be computed for example using Remez' algorithm, which is available in Sollya (http://sollya.gforge.inria.fr/); see also [13, §3.5], [START_REF] Lauter | Arrondi correct de fonctions mathématiquesfonctions univariées et bivariées, certification et automatisation[END_REF], [START_REF] Chevillard | Évaluation efficace de fonctions numériquesoutils et exemples[END_REF]. Continuing the previous example where k = 32, it turns out that a(t) = 8 i=0 a i t i can be chosen such that a 0 = 1 and a i = (-1)

i+1 A i • 2 -31 , 1 ≤ i ≤ 8, (25)
with the following values for the A i 's (including Each of them can be stored using only 32 bits and four of them are powers of two. Notice also that

A 0 = a 0 • 2 31 ): A 0 = 2 31 , A 1 = 2 30 , A 2 = 2 28 , A 3 = 134124516, A 4 =
A 8 ≤ • • • ≤ A 2 ≤ A 1 ≤ A 0 . (26)
A certified supremum norm computation (implemented for example in Sollya; see also [START_REF] Chevillard | A certified infinite norm for the implementation of elementary functions[END_REF], [START_REF] Chevillard | Certified and fast computation supremum norms of approximation errors[END_REF]) applied to this particular polynomial gives a bound less than 2 -25.5 as required by [START_REF] Fisher | Embedded Computing: A VLIW Approach to Architecture, Compilers and Tools[END_REF]. (The computed bound has the form 2 -25.972... .)

Writing an evaluation program

In our context, the evaluation program P is typically a piece of C code that implements a finite-precision computation of P (s, t). It should be accurate enough in the sense of [START_REF] De Dinechin | Assisted verification of elementary functions using Gappa[END_REF] and, since we favor latency (rather than throughput, for example), as fast as possible. Such an implementation will not require using floating-point arithmetic, and fixed-point arithmetic will suffice to evaluate P (s, t) accurately enough. In addition, we have the lemma below, which shows that P (s, t) lies in a fairly small range.

Lemma 2:

If (s, t) ∈ S × T then P (s, t) ∈ (1, 2). Proof: Let = 2 -p-3/2
. It follows from the definition of α(a) and the bound in [START_REF] Fisher | Embedded Computing: A VLIW Approach to Architecture, Compilers and Tools[END_REF] that

√ 1 + t -< a(t) < √ 1 + t + . Using 0 ≤ t ≤ 1 -2 1-p and √ 2 -2 1-p ≤ √ 2 -2 , we get 1 -< a(t) < √ 2 -. (27) Since 1 ≤ s ≤ √ 2 and < 2 -p-1 , we deduce that P (s, t) = 2 -p-1 + s • a(t) belongs to (1, 2).
With α(a) ≤ 2 -25.972... as in the previous paragraph, a sufficient condition on ρ(P) for ( 22) to be satisfied is

ρ(P) < 2 -25 - √ 2 • 2 -25.972... = 2 -26.840... .
Rounding the evaluation point. When designing an evaluation program P that achieves this accuracy, a preliminary step is to make the input (σ, τ ) machinerepresentable. On the one hand, the binary expansion of τ is 0.m λ+1 • • • m p-1 and thus, since λ is non-negative, τ is already representable using k bits provided p -1 ≤ k. On the other hand, writing RN k for rounding-to-nearest in precision k, we shall replace σ defined in (6) by

σ = 1, if e is even, RN k ( √ 2), if e is odd. (28)
The lemma below gives an upper bound on the loss of accuracy that occurs when rounding (σ, τ ) to ( σ, τ ). Lemma 3: Given p, k, s, t, a, P , P as above, let S = {1, RN k ( √ 2)} and define

ρ(P) = max (s,t)∈ S×T |P (s, t) -P(s, t)|.
Then ρ(P) < 2 1/2-k + ρ(P).

Proof: By definition, the bound ρ(P) is attained for some (s 0 , t 0 ) ∈ S × T . Writing ŝ0 = RN k (s 0 ), one has P(ŝ 0 , t 0 ) = P(s 0 , t 0 ) and, by the triangle inequality,

ρ(P) ≤ |s 0 -ŝ0 | • |a(t 0 )| + ρ(P).

The conclusion then follows from |s

0 -ŝ0 | ≤ 2 -k and the fact that (27) implies |a(t 0 )| < √ 2.
Applying Lemma 3 with k = 32 shows that we are left with finding an evaluation program P such that

ρ(P) ≤ 2 -26.840... -2 -31.5 = 2 -26.898... . (29)
Designing an evaluation program. By evaluation program, we mean a division-free straight-line program, that is, roughly, a set of instructions computing P (s, t) from s, t, p and the a i 's by using only additions, subtractions, and multiplications, without branching. In our context we shall assume first unbounded parallelism and some latencies for addition/subtraction and multiplication. Then we parenthesize the expression P (s, t) in order to expose as much ILP as we can and, thus, to decrease the overall latency of polynomial evaluation. An example of such a parenthesization, generated automatically in the same spirit as [START_REF] Harrison | The computation of transcendental functions on the IA-64 architecture[END_REF], is

P (s, t) = 2 -p-1 + s • (a 0 + a 1 t) + a 2 • (s • t 2 ) +a 3 t • (s • t 2 ) + t 2 • (s • t 2 ) • (a 4 + a 5 t) + t 2 • (a 6 + a 7 t) + a 8 t 2 .
Note that t 2 and s • t 2 are common subexpressions. With unlimited parallelism and latencies of 1 for addition, of 3 for (pipelined) multiplication, and of 1 for multiplication by a power of two (that is, a shift), such a parenthesization gives a latency of 13 (compared to 34 for Horner's scheme). We have found no parenthesization of smaller latency.

Accuracy issues. In our example, the rounding error of the program P that implements in fixed-point arithmetic the above parenthesization must satisfy (29). This requirement will be checked in Section 4.4.

For now, let us notice that this parenthesization can in fact be implemented using 32-bit unsigned integers only, which avoids to loose one bit of precision because of the need to store the sign of the coefficients a i . Indeed, an appropriate choice of arithmetic operators can be found, that ensures that all intermediate variables are positive:

P (s, t) = 2 -p-1 + s • (a 0 + a 1 t) -|a 2 | • (s • t 2 ) +a 3 t • (s • t 2 ) -t 2 • (s • t 2 ) • (|a 4 | -a 5 t) + t 2 • (|a 6 | -a 7 t) + |a 8 |t 2 . ( 30
)

Correction to ensure correct rounding

For each rounding-direction attribute • we will now obtain n = •( ) by correcting u in ( 16) and [START_REF] Chevillard | Certified and fast computation supremum norms of approximation errors[END_REF]. How to correct u depends on whether u is above or below . Thus the rounding algorithms below rely on either u ≥ or u > , which can both be implemented exactly (see Subsection 4.5).

Rounding to nearest

An algorithm producing n as in ( 9) is:

if u ≥ then n := truncate u after p -1 fraction bits else n := truncate u + 2 -p after p -1 fraction bits

If v p = 0 the above algorithm always returns the value u; this is the desired result, for in this case u is already a floating-point number and thus [START_REF] Chevillard | Certified and fast computation supremum norms of approximation errors[END_REF] implies [START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF]. If v p = 1 then u is at least 1 + 2 -p and is the midpoint between the two consecutive floating-point numbers u -2 -p and u + 2 -p : the former is returned when u > , the latter when u < , and ( 17) implies (9) in both cases; the case u = never happens because cannot be a midpoint.

Rounding down

An algorithm producing n as in ( 10) is:

if u > then n := truncate u -2 -p after p -1 fraction bits else n := truncate u after p -1 fraction bits

If v p = 1 then this algorithm always returns the floatingpoint number u-2 -p which, because of [START_REF] Chevillard | Certified and fast computation supremum norms of approximation errors[END_REF], satisfies [START_REF] Agarwal | Series approximation methods for divide and square root in the Power3 TM processor[END_REF] as required. If v p = 0 then u is already a floating-point number: if u > then u ≥ 1 + 2 1-p and the algorithm returns u -2 1-p , which is the floating-point predecessor of u, by truncating the midpoint u -2 -p ; if u ≤ , the returned value is u; in both cases, using ( 17) yields [START_REF] Agarwal | Series approximation methods for divide and square root in the Power3 TM processor[END_REF].

Rounding up

An algorithm producing n as in [START_REF] Pi | High-speed double-precision computation of reciprocal, division, square root and inverse square root[END_REF] is:

if u ≥ then n := truncate u + 2 -p after p -1 fraction bits else n := truncate u + 2 1-p after p -1 fraction bits If v p = 1 this algorithm always produces the floatingpoint number u+2 -p which, because of [START_REF] Chevillard | Certified and fast computation supremum norms of approximation errors[END_REF], satisfies [START_REF] Pi | High-speed double-precision computation of reciprocal, division, square root and inverse square root[END_REF] as required. If v p = 0 then u is already a floating-point number: if u ≥ , the algorithm returns u; if u < , it returns u + 2 1-p , which is the floating-point successor of u ≤ 2 -2 1-p ; in both cases, using ( 17) yields [START_REF] Pi | High-speed double-precision computation of reciprocal, division, square root and inverse square root[END_REF].

Summary: main steps of the computation of •( )

The box "Compute •( )" in Figures 1 and2 can be replaced by the ones in Figure 3 below. This figure recapitulates the main steps of the approach we have proposed in Sections 3.1 and 3.2 for deducing the correctlyrounded value •( ) from m and e.

IMPLEMENTATION DETAILS

The above square root method, which is summarized in Figures 123, can be implemented by operating exclusively on (unsigned) integers. We will now detail such an implementation, in C and for the binary32 format of [START_REF]IEEE standard for floating-point arithmetic[END_REF]. For this format the storage bit-width, precision, and maximum exponent are, respectively, k = 32, p = 24, e max = 127.

When writing the code we have essentially assumed unbounded parallelism and that 32-bit integer arithmetic is available. Additional assumptions on the way input and output are encoded and on which operators are available, are as follows.

Input and output encoding. The operand x and result r of sqrt are implemented as integers X, R ∈ {0, 1, . . . , 2 32 -1} whose bit strings correspond to the standard encoding of binary floating-point data (see [2, §3.4]). Our implementation of sqrt is thus a C function as in line 1 of Listing 2, which returns the value of R. Table 2 indicates some useful relationship between x and X that are a consequence of this encoding. (Of course the same would hold for r and R.) Also, the bit string of X must be interpreted as follows: its first bit X 31 gives the sign of x; the next 8 bits encode the biased exponent E of x as E = 7 i=0 X i+23 2 i ; the last 23 bits define the trailing significand field. In particular, if x is (sub)normal then

• the biased exponent E is related to the exponent e in (1) as follows:

E = e + 127 -[λ > 0], (31)
where [λ > 0] = 1 if x is subnormal, 0 otherwise; • the trailing significand field carries the bits of m in (2) as follows:

X 22 . . . X 0 = 0 . . . 01 λ bits m λ+1 . . . m 23 . (32)
Available operators. Besides the usual operators like +, -, <<, >>, &, |, ˆ, we assume a fast way to compute the following functions for A, B ∈ {0, 1 . . . , 2 32 -1}:

• max(A, B);

• the number nlz(A) of leading (that is, leftmost) zeros of the bit string of A;

• AB/2 32 , whose bit string contains the 32 most significant bits of the product AB. In our C codes, the operators corresponding to these functions will be respectively written max, nlz, and mul for readability. More precisely, with the ST231 target in mind, we shall assume that the latencies of max and nlz are of 1 cycle, while the latency of mul is of 3 cycles; furthermore, we shall assume that at most 2 instructions of type mul can be launched simultaneously at every cycle. (How to implement max, nlz, and mul in C is detailed in Appendix A.)

From Sections 4.2 to 4.5, the operand x will be a positive (sub)normal number. In this case, X 31 = 0 and since the result r is a positive normal number, R 31 = 0 as well. Therefore, it suffices to determine the 8 bits R 30 , . . . , R 23 , which give the (biased) result exponent D = 7 i=0 R i+23 2 i , and the 23 bits R 22 , . . . , R 0 , which define the trailing significand field of the result.

Handling special operands

For square root the floating-point operand x is considered special when it is ±0, +∞, less than zero, or NaN. Table 2 thus implies that x is special if and only if

X ∈ {0} ∪ [2 31 -2 23 , 2 32 ), that is, (X -1) mod 2 32 ≥ 2 31 -2 23 -1. (33)
All special operands can thus be detected using (33). The results required by the IEEE 754-2008 standard for the square root of such operands are listed in Table 3.

Note that there are essentially only two cases to consider: r is either x or qNaN. The first case occurs when x ∈ {+0, +∞, -0}, which, when x is known to be special, is a condition equivalent to

X ≤ 2 31 -2 23 or X = 2 31 . ( 34
)
If condition (34) is not satisfied then a quiet NaN is constructed by setting the bits X 30 , . . . , X 22 to 1 while leaving X 31 and X 21 , . . . , X 0 unchanged; this can be done by taking the bitwise OR of X and of the constant

2 31 -2 22 = (7F C00000) 16 ,
whose bit string consists of 1 zero followed by 9 ones followed by 22 zeros. Note that the quiet NaN thus produced keeps as much of the information of X as possible, as recommended in [2, §6.2]; in particular, the payload is preserved when quieting an sNaN, and if x is a qNaN then x is returned. Using the fact that the addition of two 32-bit unsigned integers is done modulo 2 32 and taking the hexadecimal values of the constants in (33) and (34), we finally get the C code shown in Listing 1 for handling special operands. In this code, notice that the four operations +, <=, ==, and | on X are independent of each other.

Listing 1

Code for handling special operands. Remark that (33) and (34) extend immediately to other standard formats like binaryk with k = 16, 64, 128, .... Thus, if integer arithmetic modulo 2 k is available and if the standard encoding is used for binaryk data, special values can be handled in the same way as shown here.

Computing the biased value of d and parity of e

By Property 2.1 the result r cannot be subnormal. Therefore, by applying (31) we deduce that the biased value

D = d + 127.
In order to compute D from X, we use first the expression of d in ( 8) and the relation (31) to obtain

D = (E -λ + [λ > 0] + 127)/2 . ( 35
)
Then, using (32) and the second and third rows of Table 2, we deduce that the number of leading zeros of X is λ + 8 when λ > 0, and at most 8 when λ = 0. Hence

λ = M -8, M = max(nlz(X), 8). ( 36
)
An immediate consequence of this is that

[λ > 0] = [M > 8]
. However, more instruction-level parallelism can be obtained by observing in Table 2 that

[λ > 0] = [X < 2 23 ] for x positive (sub)normal. ( 37
)
The formula (35) for the biased exponent D thus becomes

D = (E -M + [X < 2 23 ] + 135)/2 . ( 38
)
A possible C code implementing (38) is as follows:

Z = nlz(X); E = X >> 23; B = X < 0x800000; M = max(Z,8); C = B + 135; D = (E -M + C) >> 1;
Remark that in rounding-up mode Property 2.3 requires that the integer D obtained so far be further incremented by 1 when m = 2 -2 -23 and e is odd.

(

) 39
(This correction has also been illustrated in Figure 2.) Since (39) implies that x is a normal number, D must be replaced by D + 1 if and only if X ≥ 2 23 and the last 24 bits of X are 1 zero followed by 23 ones. An implementation of this update is thus:

d1 = 0x00FFFFFF; d2 = 0x007FFFFF; D = D + ((X >= 0x800000) & ((X & d1) == d2));
In fact, this treatment specific to rounding up can be avoided by exploiting the standard encoding as follows. Recall that n = •( ) is in [START_REF]IEEE standard for binary floating-point arithmetic[END_REF][START_REF]IEEE standard for floating-point arithmetic[END_REF] and has at most 23 fraction bits, and that we want the bit string 0R Hence, we implemented the computation of D -1 using the formula below, which directly follows from (38):

D -1 = (E -M + [X < 2 23 ] + 133)/2 . ( 40
)
This corresponds to the computation of variable Dm1 at line 7 of Listing 2. The only difference with the implementation of D given right after (38) occurs at line 6, where we perform B + 133 instead of B + 135.

Let us now turn to the parity of e in ( 5), which will be needed in Sections 4.3 and 4.5. Using (31), (36), and (37) we deduce that e is even if and only if the last bit of E -M + [X < 2 23 ] is equal to 1. Since the latter expression already appears in (38) or (40), an implementation follows immediately:

even = (E -M + B) & 0x1;
An alternative code, which uses only logical operators and unsigned integers, consists in taking the XOR of the last bit of E + [X < 2 23 ] and of the last bit of M :

even = ((E & 0x1) | B) ˆ(M & 0x1);
Listing 2 Square root implementation for the binary32 format, assuming a non-special operand and rounding to nearest.

1 uint32_t binary32sqrt(uint32_t X) 2 { 3 uint32_t B, C, Dm1, E, even, M, S, T, Z, P, Q, U, V;

Computing the evaluation point ( σ, τ )

In precision k = 32, rounding √ 2 to nearest gives the Q1.31 number 1 + 889516852 • 2 -31 . Thus σ in (28) is given by

σ = S • 2 -31 , (41)
with S the integer in [0, 2 32 ) such that S = 2 31 if e is even, and S = 2 31 + 889516852 = (B504F 334) 16 if e is odd. This integer S will be used in our code to encode σ, and its bit string has the form

S
τ = T • 2 -32 , (42)
where

T = 31 i=0 T i 2 i is the integer in [0, 2 32 ) such that T 31 . . . T 0 = m λ+1 . . . m 23 0 • • • 0 λ+9 . (43)
By ( 32) and (36), we see that T can be computed by shifting X left M + 1 positions. Since M is not immediately available, more ILP can be exposed by implementing this shift as in line 7 of Listing 2. Also, for the shift by M = max(nlz(X), 8) to be well defined here, the C standard [START_REF]International Organization for Standardization, Programming[END_REF] requires that 0 ≤ M < 32. One may check that this is indeed the case: since by assumption x is nonzero, the number of leading zeros of X is less than 32 and thus 8 ≤ M < 32.

Computing the approximate polynomial value v

Listing 3 below shows an implementation of the evaluation scheme (30) using 32-bit unsigned integers and the identities in [START_REF] Jeannerod | A new binary floating-point division algorithm and its software implementation on the ST231 processor[END_REF]. Multiplications by coefficients a power of two (like A 1 , A 2 , A 8 ) are implemented as simple shifts. Assuming a latency of 1 for additions, subtractions, and shifts, a latency of 3 for (pipelined) multiplications, and that at most 2 multiplications can be started simultaneously, this code can be scheduled in at most 13 cycles, as shown in Table 4. There the dashes '-' indicate that an instruction requires an additional slot because it uses an extended immediate; see Section 5.1.

The numerical quality of the code in Listing 3 has been verified using the Gappa software (see http://gappa. gforge.inria.fr/ and [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF], [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF]; see also [22, §4] for some guidelines on how to translate a C code into Gappa syntax). With this software, we first checked that all variables r 0 , . . . , r 21 are indeed integers in the range [0, 2 32 ). Then we used Gappa to compute a certified upper bound on the final rounding error; the bound produced is less than 2 -27.93 and thus less than the sufficient bound in (29). The Gappa script we wrote to perform this certification step is contained in the file binary32sqrt.gappa available at http://prunel.ccsd. cnrs.fr/ensl-00335792.

Implementing the rounding tests

There the only non-trivial part is to evaluate the expressions u ≥ (used when rounding to nearest and rounding up; see §3.2.1 and §3.2.3), and u > (used when rounding down; see § 3.2.2). It turns out that such comparisons can be implemented exactly by introducing three integers P , Q, Q which we shall define below by considering u 2 and 2 instead of u and .

Rounding down

According to the algorithm in Section 3.2.2 rounding down does not rely on the condition u ≥ but on the condition u > instead.

In order to implement the condition u > , let Q ∈ {0, 1} be such that Q = 1 if and only if equality P •2 -28 = u 2 occurs in (45), that is, if and only if P = U 2 • 2 -32 . The latter equality means that U has at least 16 trailing zeros. Since the bit string of U is [ 0 1 v 1 • • • v 24 0 0 0 0 0 0 ], this is equivalent to deciding whether v 15 = v 16 = • • • = v 24 = 0 or not. Hence the code below for computing Q :

Qprime = (V & 0x0000FFC0) == 0x0;
Property 4.2: The inequality u > holds if and only if the condition P >= Q + Qprime is true.

Proof: Since u and are non-negative, u > is equivalent to u 2 > 2 . We consider the two cases Q = 1 and Q = 0 separately. Assume first that Q = 1. Then, using the equality in (45) together with (46), we see that u > is equivalent to P > Q, that is, since P and Q are integers, to P ≥ Q + 1 = Q + Q . Assume now that Q = 0. In this case P • 2 -28 < u 2 and thus P ≥ Q implies u 2 > Q•2 -28 = 2 ; conversely, if u 2 > Q•2 -28 then, using the left inequality in (45), we find P + 1 > Q and thus P ≥ Q. Therefore, u > if and only if P ≥ Q + Q . Now, recalling that ∈ [1, 2), we deduce from (46) that Q < 2 30 . Hence Q + Q always fits in a 32-bit unsigned integer. Consequently, the mathematical condition P ≥ Q + Q is equivalent to the C condition P >= Q + Qprime.

Together with the algorithm of Section 3.2.2, the above property gives the following implementation of rounding down:

EXPERIMENTS WITH THE ST231 CORE

Combining the codes of the previous section leads immediately to a standard C implementation of sqrt for the binary32 format, with subnormal support and correct rounding for each rounding-direction attribute RN, RU, RD.

For validation purposes, each of these three versions has been compiled with Gcc under Linux (using the software given in Appendix A for emulating max, nlz, mul). This allowed for an exhaustive comparison, within a few minutes, against the square root functions of GNU C (glibc) 2 and GNU MPFR. 3 We also compiled these three versions with the Open 64-based ST200 family VLIW compiler from STMicroelectronics, in -O3 and for the ST231 core. After a review of the main features of the ST231 in Section 5.1, the performance results we have obtained in this context are presented and analysed in Section 5.2.

Some features of the ST231

The ST200 family of VLIW microprocessors originates from the joint design of the Lx by HP Labs and STMicroelectronics [START_REF] Fisher | Embedded Computing: A VLIW Approach to Architecture, Compilers and Tools[END_REF]. The ST231 is the most recently designed core of this family, and is widely used in STMicroelectronics SOCs for multimedia acceleration.

2. http://www.gnu.org/software/libc/ 3. http://www.mpfr.org/mpfr-current/ In this processor, that executes up to four integer instructions per cycle, all arithmetic instructions operate on the 64 32-bit register file and on the 8 one-bit branch register file.

Resource constraints must be observed to form proper instruction bundles containing 1 to 4 instructions: only one control instruction, one memory instruction, and up to two 32 × 32 → 32 multiplications of type mul are enabled. Other instructions can be freely used, but are limited to integer only arithmetic, without division.

Another specificity of this architecture is that any immediate form of an instruction is by default encoded to use small immediates (9-bit signed), but can be extended to use extended immediates (32-bit), at the cost of one instruction per immediate extension. For instance, up to two multiplications mul each using a 32-bit immediate can be encoded in one bundle. This makes the usage of long immediate constants such as polynomial coefficients very efficient from a memory system standpoint. On the contrary, the absence of a sophisticated data cache model (such as L2 cache) implies a quite important cost of accessing a data table (up to 130 cycles in the case of a data cache miss).

To enable the reduction of conditional branches, the architecture provides partial predication support in the form of conditional selection instructions. In the assembly line below, $q, $r, $s are 32-bit integer registers, $b is a one-bit branch register that can be defined through comparison or logical instructions: slct $s = $b, $q, $r This fragment of assembly code writes $q in $s if $b is true, $r otherwise. An efficient if-conversion algorithm based on the ψ-SSA representation is used in the Open64 compiler to generate partially-predicated code based on the slct instruction [START_REF] Bruel | If-conversion SSA framework for partially predicated VLIW architectures[END_REF].

The retargeting of the Open64 compiler technology to the ST200 family is able to generate efficient, dense, branch-free code for all the codes described in Section 4, requiring only the usage of one specific intrinsic to select the nlz instruction (number of leading zeros).

Performance results for square root on ST231

The performances obtained when compiling for ST231 with the ST200 compiler are summarized in Table 5. Versions of our codes that do not support subnormals have been implemented too, whose performances are within square brackets. Here are some observations:

• Thanks to if-conversion, the generated assembly codes are fully if-converted, straight-line programs. Thus, all latencies (numbers of cycles) and instruction numbers are independent of the input values. • Depending on •, the value N varies as expected: compared to RN, Section 4.5.1 suggests one more addition for RU, while Section 4.5.2 suggests 3 more instructions for RD (get Qprime and add it to Q). • On the contrary, the value of L turns out to be independent of •. In fact, with subnormal support, the value 23 is exactly the latency that can be expected from the codes of Section 4 when assuming unbounded parallelism. Indeed, the critical path consists of the following four sub-tasks: compute the value T in 3 cycles; compute the value V by bivariate polynomial evaluation in 13 cycles; round correctly (truncation, squaring, comparison, and selection) in 1 + 3 + 1 + 1 = 6 cycles; select the final result in 1 cycle. In other words, our implementations with subnormal support are scheduled optimally by the ST200 compiler. The same occurs for our versions without subnormals: those were obtained by replacing Z = nlz(X); M = max(Z,8); T = (X << 1) << M;

with T = X << 9;

The cost of T drops from 3 cycles to 1 cycle, and one may check that the critical path is defined by the same sub-tasks as before. Hence a new theoretical latency of 21 cycles, which again is achieved in practice. Finally, note that the overhead due to subnormal support is only of 2 cycles. • In all cases, the number of instructions per cycle (IPC) is larger than 2.6. A more precise description of bundle occupancy is given in Figure 4, assuming • = RN and subnormal support. Among the 23 bundles used, 4 contain four instructions and 12 contain three instructions. The instructions used for handling special operands are in light grey, while those corresponding to non-special operands are in black; instructions in dark grey are shared or used for selection. Note that, from the point of view of latency, handling special operands comes for free up to the last cycle, used for selecting the final result.

Table 6 shows the influence of the evaluation scheme used for P (s, t) = 2 -25 + s • a(t) on the performances of the complete square root code. Instead of evaluating P (s, t) as in (30) and Listing 3, we evaluate a(t) by three methods and then multiply by s and add 2 -25 . "Best univariate" here is one giving the lowest latency (12 cycles) we have found assuming unbounded parallelism. Comparing with the first line of Table 5, we remark that: only by 4 (less than 7%) compared to Horner's rule.

• Our scheme yields a complete code that is about twice faster than when using Horner's rule. • We use 4 cycles less than when evaluating a(t) as fast as possible and then performing in sequence the multiplication by s and the addition by 2 -25 ("Best univariate"). This shows the advantage of evaluating the bivariate polynomial P (s, t) directly, allowing for s to be distributed within a(t) = i a i t i . Table 7 gives the latencies obtained in [START_REF] Jeannerod | Faster floating-point square root for integer processors[END_REF] (for rounding to nearest and without subnormals) by implementing various other methods on ST231. We denote by R and nR the classical restoring and non-restoring algorithms; N2, G2, G1 refer to some variants of Newton-Raphson/Goldschmidt methods, based on low-degree polynomial evaluation (for the initial approximation) followed by 1 or 2 iterations; P5 and P6 are methods based exclusively on piecewise, univariate polynomial approximants (three of degree 5 or two of degree 6).

Our version in 21 cycles is about 7 times faster than the restoring method and twice faster than Goldschmidt's method G2. It is also 1.43 times faster than P6, which was the previously fastest known implementation of correctly-rounded square root on ST231.

CONCLUSION

After detailing some properties of the square root function in binary floating-point arithmetic, we proposed, for its computation with correct rounding, a high-ILP approach based on bivariate polynomial evaluation. To demonstrate the effectiveness of this approach, we further provided and analysed complete C implementations, one for each rounding-direction attribute. Some experimental results finally showed that these codes yield optimal schedules and low latencies on a VLIW integer processor like the ST231. Also, the latency overhead for subnormal support is only 2 cycles. Our approach can be extended in several ways. First, one may check that a faithfully-rounded square root follows from replacing lines 12 to 18 of Listing 2 with return (Dm1 << 23) + (V >> 7);

The latency then is 19 cycles with subnormals and 17 cycles without subnormals, leading to a speed-up by a factor of more than 1.2 in both cases.

Second, our approach is not restricted to square rooting and can be adapted to other operators: it has already been employed4 on ST231 for accelerating division [START_REF] Jeannerod | A new binary floating-point division algorithm and its software implementation on the ST231 processor[END_REF], reciprocal square root [START_REF] Jeannerod | Optimizing correctly-rounded reciprocal square roots for embedded VLIW cores[END_REF], and more generally x 1/n for |n| a "small" positive integer [START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -polynomial evaluation-based algorithms and certified code generation[END_REF]. In fact, our approach is restricted to neither algebraic functions nor the bivariate case: it could in principle be used as soon as argument reduction yields an identity of the form = F (s, t, . . .), for some multivariate function F . But it remains to determine which functions other than the ones above could benefit significantly from our approach.

Third, although we have focused in this paper on implementations for the binary32 format (single precision) and 32-bit architectures, similar codes could be derived for the binary64 format (double precision) and 64-bit architectures. The only major modification would be a new polynomial approximant together with a new bivariate evaluation scheme. However, for application reasons, we are currently mostly interested in supporting double precision on ST231. How to best use the underlying 32bit integer arithmetic available on this chip in order to emulate efficiently a binary64 square root seems to be nontrivial. It thus remains to investigate whether the approach we have introduced here can be adapted to that context.

Fig. 1 .

 1 Fig. 1. Square root algorithm for • ∈ {RN, RD}.

Fig. 2 .

 2 Fig. 2. Square root algorithm for • = RU.

  82769556, A 5 = 53306947, A 6 = 29806269, A 7 = 11452029, and A 8 = 2 21 . (See Listing 3 for their hexadecimal values). These integers A i were found by truncating the coefficients obtained after several calls to Sollya's Remez algorithm and by favoring powers of two.

x = m • 2 eFig. 3 .

 23 Fig. 3. Computation of •( ) for • ∈ {RN, RD, RU}.

  30 . . . R 23 R 22 . . . R 0 of the result r. Instead of concatenating the bit string R 30 . . . R 23 of D and the bit string R 22 . . . R 0 of the fraction field of r, one can add to (D -1) • 2 23 the integer n • 2 23 : • If n = (1.n 1 . . . n 23 ) 2 then this addition corresponds to (D -1) • 2 23 + 2 23 + (0.n 1 . . . n 23 ) 2 • 2 23 and since no carry propagation occurs, it simply concatenates the bit string of D and the bit string n 1 . . . n 23 . • If n = (10.0 . . . 0) 2 then this addition corresponds to (D -1) • 2 23 + 2 • 2 23 = (D + 1) • 2 23 . Hence R encodes the normal number r = (1.0 . . . 0) 2 • 2 d+1 .

10 V 18 return

 1018 = biv_poly_eval(S, T); // Bivariate polynomial evaluation: S [1.31], T [0.32], V [2.30] (Dm1 << 23) + ((U + 0x00000040) >> 7); 19 }

if(P >= Q) 16 return

 16 (Dm1 << 23) + ((U + 0x00000040) >> 7); else 18 return (Dm1 << 23) + ((U + 0x00000080) >> 7);

15

  Qprime = (V & 0x0000FFC0) == 0x0; 16 if(P >= Q + Qprime) 17 return (Dm1 << 23) + ((U -0x00000040) >> 7); 18 else 19 return (Dm1 << 23) + (U >> 7);

Fig. 4 .

 4 Fig. 4. Typical bundle occupancy on ST231.

TABLE 1

 1 Numerical estimations of α * d for 5 ≤ d ≤ 10.

	d	5	6	7	8	9	10
	-log 2 (α * d )	19.58 22.47 25.31 28.12 30.89 33.65

TABLE 3

 3 Square root results for special operands.

	Operand x	+0 +∞ -0 less than zero	NaN
	Result r	+0 +∞ -0	qNaN	qNaN

TABLE 2

 2 Relationship between floating-point datum x and its encoding into integer X = 31 i=0 X i 2 i . 000000000X 22 . . . X 0 with someX i = 1 [2 23 , 2 31 -2 23 ) positive normal number 0 X 30 X 29 X 28 X 27 X 26 X25 X 24 X 23 31 , 2 31 + 2 23 ) negative subnormal number 100000000X 22 . . . X 0 with some X i = 1 [2 31 + 2 23 , 2 32 -2 23 ) negative normal number 1 X 30 X 29 X 28 X 27 X 26 X 25 X 24 X 23 32 -2 23 , 2 32 -2 22 ) sNaN 1111111110X 21 . . . X 0 with some X i = 1 [2 32 -2 22 , 2 32 ) qNaN 1111111111X 21 . . . X 0

	Value or range of integer X	Floating-point datum x	Bit string X 31 . . . X 0
	0	+0	00000000000000000000000000000000
	(0, 2 23 )	positive subnormal number
				X 22 . . . X 0
			not all ones, not all zeros
	2 31 -2 23	+∞	01111111100000000000000000000000
	(2 31 -2 23 , 2 31 -2 22 )	sNaN	0111111110X 21 . . . X 0 with some X i = 1
	[2 31 -2 22 , 2 31 )	qNaN	0111111111X 21 . . . X 0
	2 31	-0	10000000000000000000000000000000
	(2 not all ones, not all zeros	X 22 . . . X 0
	2 32 -2 23	-∞	11111111100000000000000000000000
	(2

D of the exponent d of r is always given by

  Therefore, since 2 31 + 2 30 -1 = (BF F F F F F F ) 16 and given the value of the integer even (see Subsection 4.2), computing S can be done as shown at line 8 in Listing 2.The number τ in (19) satisfies τ = (0.m λ+1 . . . m 23 ) 2 . Therefore, it can be viewed as a Q0.32 number

31 S 30 S 29 . . . S 0 = 100 . . . 0, if e is even, 10 * . . . * , if e odd. Since S 31 S 30 = 10 in both cases, selecting the right bit string can be done by taking the bitwise AND of the constant (B504F 334) 16 = (10 * . . . * ) 2 and of 2 31 + 2 30 -1 + [e is even] = 110 . . . 0, if e is even, 101 . . . 1, if e is odd.

TABLE 5

 5 Performances of our approach on ST231.

			Subnormal numbers [not] supported
	•	Latency L Number N of instructions IPC = N/L
	RN	23 [21]	62 [56]	2.70 [2.67]
	RU	23 [21]	63 [57]	2.74 [2.71]
	RD	23 [21]	65 [59]	2.83 [2.81]

TABLE 6

 6 Performances with other evaluation schemes on ST231.

		• = RN and with [without] subnormals
	Evaluation scheme	L	N	IPC = N/L
	Horner's rule	44 [42] 58 [52]	1.32 [1.24]
	Estrin's method	29 [27] 60 [54]	2.07 [2.00]
	Best univariate	27 [24] 62 [56]	2.30 [2.33]

TABLE 7

 7 Latencies obtained with other square root methods on ST231, assuming • = RN and no subnormal support.

	Method	R	nR N2 G2 G1 P5 P6
	Latency L 148 133 45	42	36	33 30

The first draft of this article is[START_REF] Jeannerod | Computing floating-point square roots via bivariate polynomial evaluation[END_REF] and dates back to October 2008.

ACKNOWLEDGMENTS

This work was supported by "P ôle de compétitivité mondial" Minalogic and by the ANR project EVA-Flo. Thanks also to Nicolas Jourdan for several interesting discussions and for suggesting the use of Equation (33).

Listing 3 Bivariate polynomial evaluation code.

// A0 = 0x80000000; A1 = 0x40000000; // A2 = 0x10000000; A3 = 0x07fe93e4; // A4 = 0x04eef694; A5 = 0x032d6643; // A6 = 0x01c6cebd; A7 = 0x00aebe7d; // A8 = 0x00200000; static inline uint32_t biv_poly_eval(uint32_t S, uint32_t T) { uint32_t r0 = T >> 2; uint32_t r1 = 0x80000000 + r0; uint32_t r2 = mul(S, r1); uint32_t r3 = 0x00000020 + r2; uint32_t r4 = mul(T, T); uint32_t r5 = mul(S, r4); uint32_t r6 = r5 >> 4; uint32_t r7 = r3 -r6; uint32_t r8 = mul(T, 0x07fe93e4); uint32_t r9 = mul(r5, r8); uint32_t r10 = r7 + r9; uint32_t r11 = mul(r4, r5); uint32_t r12 = mul(T, 0x032d6643); uint32_t r13 = 0x04eef694 -r12; uint32_t r14 = mul(T, 0x00aebe7d); uint32_t r15 = 0x01c6cebd -r14; uint32_t r16 = r4 >> 11; uint32_t r17 = r15 + r16; uint32_t r18 = mul(r4, r17); uint32_t r19 = r13 + r18; uint32_t r20 = mul(r11, r19); uint32_t r21 = r10 -r20; return r21; } Truncating v = V • 2 -30 after 24 fraction bits yields

with U the integer in [0, 2 32 ) whose bit string is

Let P be the integer in [0, 2 32 ) such that

It follows from (44) and the definition of mul that

and can be represented exactly with 24 -λ bits. Since 24 -λ ≤ 32, several encodings into a 32-bit unsigned integer are possible. Because of (45) and the need to compare 2 with u 2 , a natural choice is to encode

An implementation of the computation of Q using the integer T defined in (42-43) and the parity of e can be found at line 14 of Listing 2.

Rounding to nearest and rounding up

Once the values of P and Q are available, the condition u ≥ used when • ∈ {RN, RU} can be evaluated thanks to the following characterization: Property 4.1: The inequality u ≥ holds if and only if the condition P >= Q is true.

Proof: Since u and are non-negative, u ≥ is equivalent to u 2 ≥ 2 . If u 2 ≥ 2 then, by (46) and the left inequality in (45), we deduce that P + 1 > Q, which is equivalent to P ≥ Q for P and Q integers. Conversely, if P ≥ Q then multiplying both sides by 2 -28 gives P • 2 -28 ≥ 2 and, using the right inequality in (45), u 2 ≥ 2 . To sum up, u ≥ if and only if P ≥ Q, that is, if and only if the C condition P >= Q is true.

Together with the algorithm of Section 3.2.1, this property accounts for the implementation of rounding to nearest at lines 14 to 18 in Listing 2.

Since rounding up depends on the condition u ≥ as well (see Section 3.2.3), this rounding mode can be implemented easily: simply replace lines 15 to 18 in Listing 2 with the following code fragment:

APPENDIX A SOFTWARE IMPLEMENTATION OF THE max, nlz, AND mul OPERATORS

Maximum of two unsigned integers

Number of leading zeros of a nonzero unsigned integer

Higher half of a 32-bit integer product


J Goyon

F Bertrand

O Pitois

G Ovarlez

Shear induced drainage in foamy yield-stress fluids

Keywords: 83.50.Iz, 47.56.+r, 47.50.-d

Shear induced drainage of a foamy yield stress fluid is investigated using MRI techniques.

Whereas the yield stress of the interstitial fluid stabilizes the system at rest, a fast drainage is observed when a horizontal shear is imposed. It is shown that the sheared interstitial material behaves as a viscous fluid in the direction of gravity, the effective viscosity of which is controlled by shear in transient foam films between bubbles. Results provided for several bubble sizes are not captured by the R 2 scaling classically observed for liquid flow in particulate systems, such as foams and thus constitute a remarkable demonstration of the strong coupling of drainage flow and shear induced interstitial flow.

Furthermore, foam films are found to be responsible for the unexpected arrest of drainage, thus trapping irreversibly a significant amount of interstitial liquid.

applications (transport, vibrations). From a fundamental point of view, the interdependency of drainage flow and interstitial flow induced by shear is still an open issue. Although significant progress has been done in the modelling of drainage for static foams, a more general theory, i.e. including the effects of shearing on drainage, is required. Here, we focus on the vertical (gravity) drainage behavior of a model foamy material subjected to a controlled shear in the horizontal plane. The continuous phase is a model complex liquid exhibiting a yield stress.

Whereas, due to the continuous phase yield stress, the bubble assembly is jammed at rest, we find that an imposed horizontal shear induces a fast vertical drainage. Then, we show that the coupling of shear and drainage is by far stronger than expected and that the resulting drainage behavior defies the classical results of liquid flow in particulate systems.

We study the stability of wet, monodisperse, foamy, yield-stress fluids of initial bubble volume fraction

% 1 71 ± = Φ
. Three bubble radii are studied: 125, 275 and 500 µm. The interstitial yield stress fluid is an oil in water emulsion, composed of dodecane oil droplets stabilized by surfactants described below. The droplets are 2 ± 1µm in diameter. Nearly monodisperse foams are produced by blowing a mixture of nitrogen with a small amount of perfluorohexan (C 6 F 14 ) in the surfactant solution through a porous glass frit, the porosity of which fixes the bubble size. C 6 F 14 strongly reduces the coarsening rate [10,13 ]. Two surfactant solutions are used, both in the emulsion and the foam: the first one is a solution of TetradecylTrimethylAmmonium Bromide (TTAB) at 5g/L, the second one is the same TTAB solution with dodecanol (DOH), at 0.4g/L in the foam and 2g/L in the emulsion which is high enough to saturate the oil droplets. TTAB/DOH interfaces have been shown to be less mobile (more rigid) than pure TTAB ones, thus reducing the drainage flow velocity in aqueous foams [6,14]. In the following, TTAB and TTAB/DOH interfaces will be referred to as 'mobile' and 'rigid' interfaces respectively. The foamy yield stress fluid is obtained by mixing softly the foam and the emulsion, both sharing the same continuous phase to allow their mixing without changing the bubbles and droplets interfacial properties, without any additional air introduction. The incorporation of foam into the emulsion changes the emulsion water content; then, depending on the foam water content, we may add controlled amounts of surfactant solution to the mixture to ensure that in all cases, the oil volume fraction of the interstitial emulsion is 75%, and its yield stress is 13Pa.

In order to study the stability at rest and under shear of the foamy yield stress fluids, we carry out our experiments in a coaxial cylinders Couette rheometer (inner cylinder radius: r i =4.1 cm; outer cylinder radius: r o =5 cm; height of sheared fluid: H=11 cm). The use of a wide gap allows studying large bubbles; it results here in moderate stress heterogeneity of 50%. The inner cylinder is rotating, and both cylinders are covered with sandpaper to avoid wall slip. This is the best choice that can be made for the purpose of our study as it imposes shear in the plane perpendicular to drainage (which is thus decoupled from the shear flow).

The Couette rheometer is inserted in a Magnetic Resonance Imaging (MRI) setup described in [15,16]. The bubble volume fraction is obtained both in radial and vertical directions from density imaging with an accuracy of 0.2% [17,18]. We study the time evolution of the vertical bubble concentration profiles Φ (z) for inner cylinder rotational velocity Ω ranging between 15 and 45rpm, which is high enough to ensure that in all cases, the whole gap is sheared. The local shear rate ( )

( / )/ r r V r r γ = -∂
∂ at a radius r in the gap can be deduced from the velocity profiles V(r) (Fig. 1b inset). In the following, the results are presented as a function of the macroscopic shear rate macro γ , which is the spatial average of ( ) r γ . We vary macro γ between 6 and 18s -1 .

In Fig. 1a, we first observe that the foamy yield stress fluid is perfectly stable at rest for more than 12 hours. The bubble volume fraction remains constant in time everywhere in the cell. This is consistent with the stability criterion determined for a single bubble embedded in a yield stress fluid [19]: spherical bubbles of 500µm are not expected to move as long as the yield stress of the material is higher than 2.3Pa. This result shows that the yield stress of the interstitial emulsion we use is sufficient to stabilize the foam.

Then we start shearing the material at a given macroscopic shear rate macro γ . Since the bubbles do not move in the material at rest, we have measured MRI density profiles at rest after several flow durations. The evolution of the bubble volume fraction profile of a foamy yield stress fluid under shear is shown in Fig. 1a. After introduction in the Couette cell, the sample is homogeneous. When the sample is sheared, bubbles move upwards: Φ increases up to a maximum max Φ at the top of the cell and decreases down to 0 at the bottom, the central part remaining at 0 Φ , as classically observed in creaming experiments. This bubble rise means that, under shear in the parallel plane, there is no more yield resistance to flow in the vertical direction. We checked that the bubble size and emulsion properties are not changed during the experiments: our observations are the consequence of intrinsic properties of our systems and not artifacts such as bubbles coalescence or oil droplets coalescence. A similar effect was observed for the sedimentation of solid particles in yield stress fluids [20] and granular materials [21]. After some time (of order 1000s in Fig. 1a), a steady state is finally reached: the profiles do not evolve anymore in time. The final sample is composed of pure emulsion at the bottom ( Φ =0%) and of a homogeneous dry foamy emulsion at the top ( 93 max ≈ Φ % in Fig. 1a). Moreover, we observed that max Φ depends on the bubble size and surprisingly, the amount of emulsion remaining indefinitely in the foam at this final stage can reach 10% of the total volume (Fig. 1b). Note that a careful inspection of velocity profiles in the dry part (see the inset of Fig. 1b), shows that the whole sample is sheared, so the quantity of emulsion trapped into the foam has to be attributed to an intrinsic retention mechanism.

In order to study in more details the role of the rheological properties of the interstitial yield stress fluid, it is worth studying the kinetics of drainage. From the volume fraction profiles ( , ) z t Φ , we observe a well defined front for the interstitial emulsion flow; we determine its position as p(t) where max 0

( , ) ( )/ 2 h t Φ = Φ +Φ
(see Fig. 1a). Note that we checked that the following analysis is not affected by this choice.

max 0 0 d ( ) d f V p t t = Φ -Φ Φ .
In the following, we study the impact of several control parameters on this velocity: the shear rate, the bubble size and the interface mobility.

The evolution of the front velocity V f with the macroscopic shear rate macro γ is shown in Fig. 2.

For given bubble size and surfactant solution, we observe that V f increases basically linearly with macro γ . To explain this feature, we have to understand how the flow of the interstitial fluid in gravity direction depends on shear, which is orthoradial, i.e. to analyse the 3D rheological behavior of the emulsion [20]. Emulsions behave as Herschel-Bulkley fluids [18]; a 3D form of this behavior has been proposed [20,[START_REF] Coussot | Rheometry of Pastes, Suspensions and Granular Materials[END_REF] ( ) To be more quantitative, we now compare the measured drainage velocity with that measured for liquid flow in unsheared bubbly systems. As discussed above, the interstitial liquid drained through the foam channels behaves as a purely viscous fluid, the viscosity of which depends on the shear rate in the foam channels channels γ . This latter can be estimated from the relative velocity of the bubbles, . The velocity for gravity drainage in unsheared bubble systems can be expressed as

2 n ij ij c d kd d τ τ = + , 1 2 2
( ) 2 th channels V R α η = Φ [23],
where ( ) α Φ is a coefficient depending on both the gas volume fraction and the surface mobility. Recent experiments [START_REF] Rouyer | accepted for publication in Physics of fluids[END_REF] have provided values for α for both solutions investigated in the present study: α = 17.8 and 7 kg⋅m -2 s -2 at 0.71 Φ = for mobile and rigid interfaces respectively. For macro γ = 14 s -1 , the effective viscosity of the interstitial liquid is channels η ≈ 0.15 Pa⋅s and predicted values of f V for mobiles interfaces are: 1.8 10 -3 mm/s for R = 125µm, 9 10 -3 mm/s for R = 275µm and 3 10 -2 mm/s for R = 500µm. The predicted drainage velocity compares well with the one measured for the largest bubbles (fig. 2a), but the foam made of smallest bubbles appears to drain much faster than expected, at least by one order of magnitude. Surprisingly, we actually observe in Fig. 2a that the bubble radius has very limited impact on the drainage velocity f V . In aqueous foams, f V scales as R 2 [14] and would change by a factor 16 from the smallest to the largest bubbles; here, it increases by only a factor 1.5.

From now, we propose to interpret these results by considering the flow behavior of the yield stress fluid in foam films. It has been shown that drainage of an isolated horizontal foam film containing colloidal suspensions [START_REF] Sethumadhavan | [END_REF] leads to a decrease of the film thickness down to a critical value below which at least one layer of particles is trapped in the film. In the sheared foamy emulsion, transient films are formed from collisions of neighbouring bubbles [25] and we expect the behavior of those films to be similar to that reported for an isolated film [START_REF] Sethumadhavan | [END_REF].

Thus, it is expected that in the foamy emulsion, foam films always contain a few trapped layers of emulsion droplets that cannot be ejected, thus defining a minimum film thickness f h .

The minimum volume fraction associated to those foam films at the end of drainage can be determined using the geometrical model of the well-known Kelvin cell: the films' surface area per bubble is 2 14R S f ≈ [7] and the corresponding film volume is 2 / 2 7

f f f f V S h h R = ≈
. The maximum gas volume fraction thus writes max 1 1/(1 0.6 / )

f R h Φ ≈ - +
. As shown in Fig. 1b, a very good agreement is found with experimental values assuming f h = 10µm, which indeed corresponds to a few emulsion droplet diameters. Note that this large value differs by one or two orders of magnitude from the one measured for aqueous foams in draining conditions [9].

We argue that this strikingly high value for the film thickness is also responsible for the surprising drainage behavior reported above. Under shear, during the approach of two neighbouring bubbles, the interstitial liquid is squeezed out from the contact area so that the film thickness decreases at any stage of its existence from its initial value h 0 down to the final minimum value h f . During the film thinning process, a significant amount of emulsion is then drained towards connected foam channels, thus increasing the effective flow rate inside those channels. To go one step further, one has to evaluate the emulsion effective viscosity in the films, i.e. to estimate the shear rate films γ in the films. For dry aqueous foams, shear results in both a squeeze flow and a shear flow in the films [25]; the squeeze flow results in an important dynamic pressure in the films that has to be counterbalanced on the film surface by the bubble capillary pressure, which thus controls (and limits) the squeeze velocity (see Eq. 1 in [25]); it can then be shown that dissipation occurs mainly in the shear flow [25]. Here this is not the case. It should be noted that our systems are initially very wet and that, as shown above, the film are thick (h>10µm). Then, the capillary pressure 2 / R γ (from 640 Pa for the 125µm bubbles to 160 Pa for the 500µm, γ is the surface tension) has to be compared to the pressure due to the squeeze flow of the interstitial yield stress fluid (the emulsion) at yield [START_REF] Coussot | Rheometry of Pastes, Suspensions and Granular Materials[END_REF]:

2 3 f c R h τ , where 1 2 ( ) f R Rh ≈
is the typical size of the film area; this dynamic pressure varies from 30 Pa for the 125µm bubbles to 60 Pa for the 500µm, if h=h f =10µm. This shows that the velocity of the squeeze flow in our systems is not limited by the capillary pressure and thus that the decrease in thickness of a fresh foam film is initially controlled by the relative velocity of the two bubbles, i.e. . The average liquid velocity at the film/channel transitional area is thus given by: ( )

2 macro f u R R h γ ≈
, where h an average thickness value during the film thinning, expected to be slightly above h f , i.e. for the 125µm bubbles for example.

Obviously, the coupling of flows in films and channels is efficient only if the volume of emulsion drained from one film is comparable to that contained in one channel. The film to channel volume ratio is ( )( )

3 2 2 2 2 f m R h R h R δ δ ≈ ≈
(where ≈ . This is in good agreement with measured velocities. As a summary, while shearing has no significant effect on the dynamics of aqueous foam drainage [27], it drastically modifies the drainage behavior of the complex liquid: there is an additional contribution of the films on drainage, which increases with decreasing bubble size.

3 2 δ π = - [ 26
Finally, the effect of surfactant mobility is considered. A comparison of results obtained for 'mobile' and 'rigid' interfaces is presented in Fig. 2b. Surprisingly, the impact of the surfactant on the drainage velocity is opposite to what is usually measured for 'static' aqueous foams [14,28]. Indeed, we observe that drainage is more than 2 times faster for 'rigid' interfaces than for 'mobile' ones, for a given macroscopic shear rate. Qualitatively, this effect can be understood by considering that the shear between bubbles (films and channels) is more efficient for 'rigid' interfaces than for 'mobile' interfaces, leading to a further decrease of the interstitial fluid effective viscosity. Additionally, the surface tension of TTAB/DOH solution is lower than for pure TTAB [10], meaning that bubbles are more easily flattened during collisions and that the discussed above coupling of flows is stronger for TTAB/DOH.

In conclusion, we have studied the stability under gravity of model wet foamy yield stress fluids. We have shown that while the yield stress of the interstitial fluid stabilizes the foam at rest, rapid drainage is induced by shear. The interstitial yield stress fluid then behaves as a viscous fluid in the direction orthogonal to shear, which is that of gravity, and is thus drained.

Striking differences with the case of static aqueous foams have been observed, namely the independence of the velocity on the bubble size and the unexpected arrest of drainage at high liquid fraction: these last behaviors can be understood by considering the shear induced flow of interstitial material in the transient foam films and its coupling with the one in the foam channels.

The authors thank Institut Carnot VITRES and the French Space Agency for funding.

Fig. 1 a

 1 Fig. 1 a) Bubble volume fraction Φ vs. height z of a foamy yield stress fluid with 275µm radius bubbles, for a 12 hours rest (*), and as function of shear duration at macro γ

  shear rate, ij d is the strain rate tensor and c τ the yield stress of the material. From this expression, it can be seen that the shearing in any direction affects the shear stresses in all directions. If the shear rate in the interstitial emulsion due to shear in the parallel plane, but dominates over the strain rate due to drainage, then the shear resistance rz τ to a small vertical strain rate rz d induced by drainage reads with our observation of a drainage velocity proportional to macro γ . Consistently, it has been recently shown that the sedimentation velocity of a dilute suspension of hard spheres in a sheared yield stress fluid scales with the Stokes velocity in a Newtonian fluid of viscosity eff η [20]. Finally, this analysis allows the comparison of drainage velocities with expressions valid for interstitial Newtonian media, the relevant viscosity to consider problem in the following is to properly evaluate interstitial γ in connection with the foam geometry.

Fig. 2 :

 2 Fig. 2: Drainage front velocities vs. macroscopic shear rate. (a) Influence of the bubbles radius: (■) 125µm, (•) 275µm and (▲) 500µm. (b) Influence of the surface mobility for the 275µm bubbles: (•) TTAB and (♦) TTAB+DOH.

  so that only a few percent of the liquid in the channels comes from the films: this explains why the measured drainage velocities are in agreement with the first theoretical analysis based on channels γ. In opposition, for the 125µm bubbles, nearly 50% of the liquid contained in the channels comes from the films. As a result, the apparent viscosity of the liquid in the channels is mainly controlled by films γ , and the previous theoretical drainage velocity based on channels γ

  By analogy with classical drainage experiments, we focus on the drainage velocity:

	(	)(	)

The bottom of the Couette cell can have disruptive effects, so we analyze the drainage kinetics only in its first stage (when p evolves by less than 15% of the sample's total height). In this case, p(t) varies linearly in time, as classically observed in creaming experiments, which allows us to define properly a front velocity.