A Al-Zein

Ph Papet

O Mathon

B Hehlen

C Levelut

J Haines

J Rouquette

Role of Zn and Nb in the Giant Piezoelectric Response of PbZn 1/3 Nb 2/3 O 3

Relaxor ferroelectrics perovskite are highly polarizable and can exhibit giant coupling between elastic strain and an applied electric field. Here, we report an "in situ" extended X-ray absorption fine structure (EXAFS) study of a PbZn 1/3 Nb 2/3 O 3 (PZN) single crystal as a function of the electric field. We show that the strong dipoles in the NbO 6 octahedra bonds are aligned along the 4 <011> directions close to the orientation of the electric field, while a small reversible polar shift occurs for Zn in the direction of the electric field, i.e. positive or negative. This reversible Zn-O polar shift is proposed to play an important role in both the "easy" switching of the ferroelectric polarisation and of the giant piezoelectric effect in PZN.

Introduction

The growth of PbZn 1/3 Nb 2/3 O 3 (PZN) and PZN-based perovskite single crystals has allowed a significant advance in high piezoelectric efficiency to be obtained, 1 which can find a large number of technological applications (actuators, sensors, microelectromechanical systems) 2 . These materials exhibit ultra-high piezoelectric coefficients d 33 from 1100 to 2500 pC/N in PZN and PZN-8%PT (a PZN solid solution with 8% of PbTiO 3 ) respectively for <001> oriented crystals 3, 4 , which can be compared with the 220 pC/N in the wellknown lead zirconate titanate ceramics (PZT and 550 pC/N for doped PZT) 5 . The origin of the large dielectric susceptibility and piezoelectric properties in these highly polarizable ferroelectric perovskites of general formula ABO 3 is due to the cooperative displacement of the B cations, i.e. the Zn/Nb atoms in our case, from the center of symmetry of their octahedral sites. These spontaneous displacements create local dipoles that are aligned parallel to one another in order to give a macroscopic dipole moment by the application of an electric field E. In the case of PbMg 1/3 Nb 2/3 O 3 and PZN single crystals, the average ABO 3 pseudocubic perovskite structure is perturbed by translational disorder on an intermediate length scale 6, 7 . It is known that this nanoscopic inhomogeneity is responsible for relaxor properties of these materials, such as smearing and frequency dependence of the dielectric anomaly associated with the ferroelectric ordering 8 and for phenomena implying the presence of polar nanoregions (PNR) persisting hundreds of degrees above the temperature of the dielectric permittivity maximum 9, 10 . While focusing on the PZN-PT solid solution, one can remark the large difference in piezoelectric properties of the two endmembers PbTiO 3 (d 33 = 35 pC/N) 11 and PbZn 1/3 Nb 2/3 O 3 (d 33 = 1100 pC/N) 1 , which have tetragonal (P4mm space group) 12 and pseudocubic structures (with small rhombohedral distortion) 13 respectively. As long range order in ferroelectric relaxor is not properly defined, it appears more appropriate to differentiate these two crystals based on their individual atomic arrangements [14][15][16][17][18] . This is why we decided to perform an "in situ" extended X-ray absorption fine structure study of a PbZn 1/3 Nb 2/3 O 3 single crystal as a function of the applied electric field to clarify the individual role of the displacements of the Zn and Nb atoms. In the PZN-PT solid solution, the role of Pb atoms is of particular importance as it hybridizes with oxygen states leading to a large strain linked to the piezoelectric effect 19 . Additionally, if lead atoms are replaced by barium atoms, i.e. BaZn 1/3 Nb 2/3 O 3 , 1:2 ordered structures (1 Zn/ 2 Nb) with a trigonal distortion space group P3 m1 are observed, which is a clear indication of the specific role of lead in the disordering on B-site positions 20 . Hence, the clear difference in the ABO 3 perovskite between PbTiO 3 and PbZn 1/3 Nb 2/3 O 3 is based on the B-site elements. In PbTiO 3 , small Ti 4+ ions (r(Ti 4+ ) = 0.605 Å) 21 exhibit pronounced off-center shifts (0.3 Å) 22 in their oxygen octahedra, giving rise to the strong spontaneous polarization (P S = 75 μC/cm 2 ) [START_REF] Gavrilya | Soviet Physics Solid State[END_REF] along the [001] tetragonal direction. Looking at the unoccupied or fully occupied d orbitals, (d 0 state of the Nb 5+ and the d 10 state of the Zn 2+ ions), it is possible to have an insight into the local shifts on the same site, i.e. the center of the oxygen octahedron. Zn 2+ and Nb 5+ have ionic radii of 0.74 and 0.64 Å 21 , respectively. One can therefore expect a much more pronounced off-center shift of smaller Nb atoms in their oxygen octahedral cage in comparison to the larger Zn atoms. This scenario was confirmed based on X-ray absorption spectroscopy (XAS) experiments; the Nb off-center displacement is experimentally found to be about 0.1 Å, while Zn exhibits no shift 24, 25 . Based on the information on B-site elements described above, one can predict that the large Tioff center shift in PbTiO 3 will require a large coercive field (6.75 kV/cm) [START_REF] Gavrilya | Soviet Physics Solid State[END_REF] to flip the spontaneous polarization, which can be linked with the relatively low piezoelectric coefficient of lead titanate. In contrast, PZN exhibits a smaller average off-center shift, which can explain the low coercive field (2.7 kV/cm) needed to reverse the smaller P S (P S = 25 μC/cm 2 ) 1 .

Hence, the B-site polar character in PZN arises predominantly from the Nb-atoms. However, the mechanism underlying the ultrahigh performance of PZN and more particularly the electric field dependence of the Zn and Nb atoms is the key issue to give the appropriate atomic scale picture of this functional material under "in-situ" conditions corresponding to those experiments in technological applications. Therefore, XAS appears as the technique of choice to clarify the individual role of the displacements of the Zn and Nb atoms as a function of the applied electric field. Here we report an "in situ" EXAFS study on PZN single crystal poled along the [001] direction, i.e. the direction of the giant piezoelectric effect.

Experimental Section

EXAFS measurements were performed on PZN plates metalized with gold (typically 2. . The contribution of interest to the total EXAFS signal is given by the single scattering between the photoabsorber, i.e. the B-atom (Zn or Nb), and the first six oxygen neighbours over the whole electric field range, i.e. 0 ≤ E ≤ 8 kV/cm. The time to measure one EXAFS dataset as a function of E was about 75 minutes. At least two spectra per electric field value were needed to give reproducible data. Attempts were carried out to measure EXAFS spectra during the first polarization cycle of the PZN crystal. However, kinetic effects during polarization are well-known in ferroelectric PZN crystals 28 , i.e. the electric field necessary to pole the sample strongly decreases with time. This means that the poling state of an unpoled sample under electric field can vary during acquisition and the set up used in this study does not allow us to accurately describe the EXAFS spectra of the first polarization cycle in PZN. Thus, we will only present measurements on the unpoled state and the poled state under electric field.

Results

High quality Nb K edge and Zn K edge EXAFS data of PZN were obtained in the 0 ≤ E ≤ 8 kV/cm electric field range, Figure 1. k 3 χ(k) spectra are shown here to better take into account the contribution of the Nb(Zn)-O first coordination shell, which is linked to the fast damping of the scattering from the O-atoms with increasing k. In the case of the Nb-EXAFS data, Figure 1a, one can observe that the electric field has no drastic effect on the positions of the EXAFS fine structure, which indicates that the Nb local environment is not strongly affected during poling. However, the oscillations exhibit stronger amplitudes at high k between 8 and 11. Let us now discuss the Zn-EXAFS data, Figure 1b. As in the spectra at the Nb K edge, one can observe that the electric field has no strong effect on the positions of the EXAFS fine structure, which is consistent with a similar local environment of Zn atoms during poling. A careful inspection of the EXAFS oscillations shows that the unpoled crystal and the poled crystal at 2.5 kV/cm exhibit higher amplitudes between 6 and 10 Å -1 compared with the poled samples at 4 and 6 kV/cm respectively. However, when the electric field is removed (poled state 0 kV/cm) the pronounced high k oscillations reappear. Again, the main effect of the electric field on the Zn K edge EXAFS signal will be based on the orderdisorder around the Zn-atoms and its associated σ All the different structural models were tested for the Nb-O first coordination shells, i.e. orthorhombic (tetragonal), rhombohedral and cubic local symmetries. We have to keep in mind that ferroelectric material can exhibit either no Bsite off-center, with the corresponding cubic local structure, or a polar distortion, rhombohedral or tetragonal; this latter tetragonal distortion can only be orthorhombic (4-2) with a corresponding dipolar moment along the <110> direction, Figure 3, as the Jahn-Teller tetragonal distortion presents no polar shift. )) which was constrained to be the same for R 1 and R 2 , and the number of neighbours, which was fixed to six. σ 2 (Nb-O) contains two contributions: a thermal (σ 2 dyn. ) and a structural (σ 2 stat. ) contribution linked to the local vibrational dynamics and the local structural distortions, respectively. In this study, as the thermal contribution is constant, changes in σ 2 (Nb-O) correspond to a structural contribution, linked with the orthorhombic distortion constrained by R 1 and R 2 .This orthorhombic Nb-O environment, Figure 3, was already found to be a common feature in all the PbB'NbO 3 perovskite compounds (B'=Zn, Mg, Sc, In) and is independent of the macroscopic symmetry (cubic, rhombohedral, orthorhombic) , which is also in agreement with our data for PbZr 0.52 Ti 0.48 O 3 30 . Figure 2 shows the backtransformed EXAFS signal corresponding to the Nb-O atom pair distribution function, i.e. backward Fourier transform, in the k 3 -weighted space calculation obtained for an unpoled sample, a sample poled at 4 kV/cm and a poled sample where the electric field has been removed (0 kV/cm). As already observed in the Nb K edge spectra 21 , one can clearly observe that the oscillations of the Nb-O atom pair distribution function between 8 and 11.3 Å -1 exhibit a higher amplitude (and identical positions) for the poled sample at 4 kV/cm and 0 kV/cm respectively, compared to the unpoled state, which is a clear indication of a change in the order in the Nb-O local environment. Figure 5a and 5b show the inverse Fourier transformed EXAFS signal along with the best fit calculation obtained for the poled sample at 4 kV/cm and 0 kV/cm after field removal respectively. The experimental data are correctly reproduced over the electric field range investigated.

Zn-O disordered model

Let us now discuss the Zn-EXAFS data. Here, we refined the k space spectra isolating the Zn-O first coordination shell between 1.2 Å and 2 Å and the back transforming in k space (2 ≤ k ≤ 10 Å -1 ) using a standard procedure. Over the entire electric field range, a symmetric distribution of bond dis-tances was used for Zn-O as no "clear distortion" was observed up to the highest electric field value as explained below. Additionally as already mentioned in the introduction, the large Zn atoms are supposed to exhibit no polar shift and this symmetric distribution of bond distances is clearly appropriate as in the case of the Zr atoms in PZT 31 . We attempted to use rhombohedral and tetragonal models to fit our data, i.e. with 2 distances R 1 and R 2 . Below 4 kV/cm R 1 =R 2 , the data are consistent with the proposed symmetric model, whereas a distortion could be observed for E ≥ 4 kV/cm although the estimated standard deviation was found to be of the order of the proposed distortion. Additionally, as the momentum transfer is limited to 10 Å -1 in our data, both rhombohedral and tetragonal structural models gave similar agreement factors. Therefore, we decided to keep the symmetric distribution of bond distances for the entire electric field range as all the data are analyzed with the same structural model. In such a disordered model as the thermal contribution is constant, the increase in σ 2 (Zn-O) thus corresponds to a polar structural contribution µZn, i.e. ΔR between R 1 and R 2 is therefore contained in the DWF. Additionally such a disordered model only uses two refined parameters R and σ 2 (Zn-O) instead of at least three for an ordered model (R 1 , R 2 , σ 2 (Zn-O)), which is particularly interesting using the restrained momentum transfer range obtained in our data. In the entire electric field range we used the two refined following structural parameters: R (average bond distance, R = 2.066 (5) Å), σ 2 (Zn-O), and the number of neighbours which was fixed to six. The backtransformed EXAFS signal for an unpoled sample, a sample poled at 6 kV/cm and a poled sample after removal of the electric field (0 kV/cm) is shown in Figure 6. 3 -weighted space calculation obtained for an unpoled sample (blue), a sample poled at 6 kV/cm (black) and a poled sample (red) after removal of the electric field (0 kV/cm).

As already observed in Zn-EXAFS data, Figure 1b, it is possible to note that the oscillations of the Zn-O atom pair distribution function exhibit a stronger amplitude (and identical positions) between 6 and 10 Å -1 for the unpoled and poled (0 kV/cm) sample, compared to the poled (6kV/cm) state, which shows a change in the order in the Zn-environment. Figures 7a and7b show the backtransformed EXAFS signal along with the best fit calculation obtained for the poled state at 0 and 6 kV/cm respectively. The experimental data are correctly reproduced over the electric field range investigated. -O) with electric field corresponds to a structural contribution (i.e. polar displacement), which cannot be detected within the experimental accuracy of the EXAFS measurements with the used limited momentum transfer range, but can clearly be seen in the increase in the static contribution, which is significantly greater than the esd. As already discussed, adding a distortion with 2 distances R 1 and R 2 could be useful above 4 kV/cm although the estimated standard deviation was found to be of the order of the proposed distortion due to the limited k-range of 2 -10 Å -1 . This indicates that this increase in σ 2 (Zn-O) as a function of the electric field using a symmetric distribution of bond distances can be linked with a small, electric-field-induced polar shift. For E ≥ 4 kV/cm, this polar structural distortion reaches its critical limit which could be defined as the saturation of this polar displacement. Additionally, this small polar shift is found to be reversible, Figure 8b. Based on the results obtained by EXAFS at the Nb and Zn K edges, we can propose a simplified scheme of the B-site behaviour under electric field in the PZN poled sample. As the electric field dependence of Zn-O and Nb-O bond lengths are found to be within the estimated standard deviation, the critical parameter appears to be the DWF, which in this study accounts for the static distribution of the corresponding Zn-O and Nb-O bonds. In the AB'B''O 3 structure of PZN, small Nb atoms are at the origin of strong dipoles along the <011> pc (with <011> pseudocubic orientations) directions. When the electric field is applied along <001> pc direction of a poled sample, Figure 9, the electric field will induce a small polar shift to the larger Zn atoms which will create a favourable energetic situation for the alignment of the strong dipoles bonds along <011> pc close to the orientation of the electric field. As this small polar shift can be obtained with a modest electric field, the corresponding E C to align the µNb dipoles is relatively low, which may account for the large piezoelectric effect in PZN. Additionally, this small polar shift of Zn is reversible, which means that it can change its direction when the applied electric field is reversed; such a reversible easy switching can clearly explain the flipping of the spontaneous ferroelectric µNb polarisation, Figure 9. In this simplified scheme of the B-site behaviour, the Zn-polar shift could be qualitatively viewed as the difference between the spontaneous polarisation at saturation and the remnant polarisation: P S(saturation) -P S(remnant) .

Conclusions

To conclude an "in situ" EXAFS study as a function of the applied electric field was undertaken on giant piezoelectric PZN single crystals at both the Nb and Zn K edges to understand the atomic scale behavior of the B-site atoms under the conditions of use in technological applications. We find that the small Nb atoms are at the origin of highly polarized Nb-O bonds, which are aligned in the poled state.

Alignment of the dipoles can be induced by a small polar shift of the Zn atoms under electric field, which plays an important role in the small E C of PZN and thus explains the giant piezoelectric response. As this small polar shift is found to be reversible as a function of the electric field, it therefore plays a major role in the polarization flipping when the direction of the electric field is reversed. Finally, such a mechanism could also be generalized in other ferroelectric materials which also present high or giant piezoelectric effect, the PbMg 1/3 Nb 2/3 O 3 (PMN), (Pb 1-y La y )Zr 1-x Ti x O 3 (P(L)ZT), BiScO 3 -PbTiO 3 solid solutions 35 ...

26 ,

 26 5 mm  2.5 mm  200 μm), with large surface perpendicular to the [001] direction, at beamline BM29 of the European Synchrotron Radiation Facility. Direct voltage was applied between the electrodes along the [001] direction during the measurements. A (311) Si monochromator with a resolution of 0.5 eV was used. Because of the large absorption coefficient of lead, EXAFS data were collected in fluorescence mode at both the Zn K edge (9659 eV) and the Nb K edge (18986 eV). The fluorescence mode appears to be an appropriate geometry to investigate the electric field dependence of the crystal, i.e. 45° scattering geometry between the X-ray incident direction electric field vector 𝑒 and the applied electric field 𝐸 linked to the polar shift. The data have been corrected to the self-absorption approximation, which affects the amplitude of EXAFS spectra measured in fluorescence mode. To avoid any correction errors, we checked that the EXAFS spectrum of a pellet of finely ground powder of PZN and boron nitride measured in transmission geometry gives similar results as the fluorescence spectrum of our unpoled PZN single crystal. It is worth noting that no changes are expected in the XANES and the pre-edge regions in this system especially the p-d peak, which arises from mixing of unoccupied p and d states of the B atom because the d orbitals of the Zn atom are fully occupied (d 10 state) and Nb K hole width is large (γK(Nb) ≈ 4.14 eV). The EXAFS data analysis was carried out with the IFEFFIT package27

Figure 1 .

 1 Figure 1. EXAFS spectra of PbZn 1/3 Nb 2/3 O 3 crystal as a function of electric field. (a) Nb K edge spectra, (b) Zn K edge spectra. NP and P accounts for non poled and poled state.

2 )

 2 samples. This can clearly be associated with the EXAFS Debye-Waller factor: exp (2σ 2 k ), which accounts for the relative mean square displacement linked to the Nb-O distances and is therefore an indication of order/disorder.

2 ,

 2 but in a different manner when compared with the Nb-data.Nb-O ordered modelTypical Nb-O first coordination shell in the k space for different electric field values are shown on Figure2. Fitting of the EXAFS signal was done in k space isolating the Nb-O first coordination shell between 1.1 Å and 2.1 Å and back transforming in k space (2 ≤ k ≤ 11.3 Å -1 ) using a standard procedure.

Figure 2 .

 2 Figure 2. Inverse Fourier transform of the Nb-O atom pair distribution function in the k3 -weighted space calculation obtained for an unpoled sample (black), a sample poled at 4 kV/cm (red) and a poled sample (blue) after the electric field was removed (0 kV/cm).

Figure 3 .

 3 Figure 3. Orthorhombic 4 short (red, R 1 = 1.94 (1) Å) -2 long (yellow, R 2 = 2.15 (1) Å) distribution of Nb-O bond distances along with <110> Nb-off center μ Nb .

Figure 4 . 2 (

 42 Figure 4. Inverse Fourier transform EXAFS signal corresponding to the Nb-O atom pair distribution function along with the best fit calculation obtained for the unpoled crystal using: (a) Orthorhombic (R 1 = 1.94 (1) Å), R 2 = 2.15 (1) Å) distribution of Nb-O bond distances. (R-factor = 0.5%), (b) rhombohedral 3 -3 (R 1 = 1.97 (4) Å, R 2 = 2.19 (4) Å) distribution of Nb-O bond distances (R-factor = 14.5 %), (c) cubic -6 (R = 1.92 (1) Å) distribution of Nb-O bond distances (R-factor = 3.6 %).This orthorhombic distribution of bond distances was therefore used, characterized by the 3 following refined structural parameters: R 1 (four equivalent short distances, R 1 = 1.94 (1) Å) and R 2 (two equivalent long bond distances R 2 =

Figure 5 .

 5 Figure 5. Inverse Fourier transform EXAFS signal corresponding to the Nb-O atom pair distribution function along with the best fit calculation obtained for (a) the poled sample at 4 kV/cm and (b) the poled sample after removal of the electric field (0 kV/cm).

Figure 6 .

 6 Figure 6. Backward Fourier transform in the k3 -weighted space calculation obtained for an unpoled sample (blue), a sample poled at 6 kV/cm (black) and a poled sample (red) after removal of the electric field (0 kV/cm).

Figure 7 . 2 (

 72 Figure 7. Inverse Fourier transform EXAFS signal corresponding to the Zn-O atom pair distribution function along with the best fit calculation obtained for (a) the poled crystal (0 kV/cm) after field removal and (b) the poled state (6 kV/cm).

Figure 8 . 2 ( 2 (σ 2 ( 2 ( 2 ( 2 (

 8222222 Figure 8. (a) Electric field dependence of the extracted σ 2 (Nb-O). Disordered (12-O) and ordered polarized dipoles (4-O) are shown in the insets 32 . (b) Electric field dependence of σ 2 (Zn-O). Filled squares correspond to points obtained on increasing field and the open circle corresponds to the point obtained after the field removed.

Figure 2 ( 2 ( 2 ( 2 )

 2222 Figure 8b shows the electric field dependence of the experimental σ 2 (Zn-O) using a disordered cubic model. A noticeable increase of σ 2 (Zn-O) can be observed with increasing E followed by a flat slope about 0.0131 Å 2 above 3 kV/cm. It is interesting to note that electric field value is close to that of the coercive field (E C ). Additionally when the electric field is removed, σ 2 (Zn-O) comes back to its initial value (0.0105 Å 2 ) while the refined Zn-O average distance remains almost constant, i.e. R = 2.066 (5) Å. In this study, as the thermal contribution is constant, changes in σ 2 (Zn-O) with electric field corresponds to a structural contribution (i.e. polar displacement), which cannot be detected within the experimental accuracy of the EXAFS measurements with the used limited momentum transfer range, but can clearly be seen in the increase in the static contribution, which is significantly greater than the esd. As already discussed, adding a distortion with 2 distances R 1 and R 2 could be useful above 4 kV/cm although the estimated standard deviation was found to be of the order of the proposed distortion due to the limited k-range of 2 -10 Å -1 . This indicates that this increase in σ

Figure 9 .

 9 Figure 9. Simplified scheme of the B-site behaviour under electric field in the PZN poled sample.

Acknowledgments

The authors would like to thank the beamline staff at BM29 at the ESRF for their assistance and Alexandre Lebon who synthetized the single crystals of PZN during its PhD thesis.

ASSOCIATED CONTENT Supporting Information

AUTHOR INFORMATION

Corresponding Author

Jerome.Rouquette@univ-montp2.fr

Author Contributions

(Word Style "Section_Content").

‡These authors contributed equally. (match statement to author names with a symbol, if applicable)

Notes

The authors declare no competing financial interests. Any additional relevant notes should be placed here.


B Zhang

M J Fadili

J-L Starck

MULTI-SCALE VARIANCE STABILIZING TRANSFORM FOR MULTI-DIMENSIONAL POISSON COUNT IMAGE DENOISING

We propose in this paper a Multi-Scale Variance Stabilizing Transform (MSVST) for approximately Gaussianizing and stabilizing the variance of a sequence of independent Poisson random variables (RVs) filtered by a low-pass linear filter. This approach is shown to be fast, very well adapted to extremely low-count situations and easily applicable to any dimensional data. It is shown that the RV transformed using Anscombe VST can be reasonably considered as stabilized for an intensity λ 10, using Fisz VST for λ 1 and using our VST (after low-pass filtering) for λ 0.1. We then use the MSVST technique to stabilize the detail coefficients of the Isotropic Undecimated Wavelet Transform (IUWT) of multi-dimensional Poisson count data. We use a hypothesis testing framework in the wavelet domain to denoise the Gaussianized and stabilized coefficients, and then apply the inverse MSVST-IUWT to get the estimated intensity image underlying the Poisson data. Finally, potential applicability of our approach is illustrated on an astronomical example where isotropic structures must be recovered.

INTRODUCTION

The ability to restore the underlying intensity from an inhomogeneous Poisson process is crucial for many applications. We observe a discrete dataset of counts x = (X n ) n∈I where I is the index set. Each count X n is independently Poisson distributed with a mean λ n , i.e. X n ∼ P(λ n ).

A host of estimation methods have been proposed in the literature. A common solution is to use a variance stabilizing transform (VST), which "Gaussianizes" the Poisson noise before applying the standard wavelet thresholding denoising on the transformed signal. For example, Anscombe transform [START_REF] Anscombe | The Transformation of Poisson, Binomial and Negative-Binomial Data[END_REF] and Fisz transform are respectively proposed in [START_REF] Donoho | Nonlinear wavelet methods for recovery of signals, densities and spectra from indirect and noisy data[END_REF] and [START_REF] Fryźlewicz | A Haar-Fisz algorithm for Poisson intensity estimation[END_REF]. Besides nonlinear VST, direct wavelet filtering has been studied in [START_REF] Nowak | Wavelet-Domain Filtering for Photon Imaging Systems[END_REF][START_REF] Antoniadis | Wavelet shrinkage for natural exponential families with quadratic variance functions[END_REF]. The state-of-the-art methods are Bayesian approaches (see overview in [START_REF] Besbeas | A Comparative Simulation Study of Wavelet Shrinkage Estimators for Poisson Counts[END_REF]). A great part of the above methods are based on the Haar wavelet transform until recently, Jansen [START_REF] Jansen | Multiscale Poisson data smoothing[END_REF] introduced a conditional variance stabilization, which generalizes the idea in [START_REF] Fryźlewicz | A Haar-Fisz algorithm for Poisson intensity estimation[END_REF] and is applicable to any family of wavelet transforms. A Bayesian scheme was also derived within this framework, which can be deemed as an extension of [START_REF] Kolaczyk | Bayesian multiscale models for Poisson processes[END_REF][START_REF] Timmermann | Multiscale Modeling and Estimation of Poisson Processes with Application to Photon-Limited Imaging[END_REF][START_REF] Nowak | A statistical multiscale framework for Poisson inverse problems[END_REF]. Other technique such as penalized maximum likelihood estimation [START_REF] Willett | Multiscale likelihood analysis and image reconstruction[END_REF][START_REF] Sardy | Automatic smoothing with wavelets for a wide class of distributions[END_REF] can also be considered Bayesian, since the penalization term implicitly introduces a prior on the underlying intensity. In general, Bayesian methods outperform those of direct wavelet filtering [START_REF] Nowak | Wavelet-Domain Filtering for Photon Imaging Systems[END_REF][START_REF] Antoniadis | Wavelet shrinkage for natural exponential families with quadratic variance functions[END_REF]. However Bayesian approaches require having enough "useful signals" in the observations in order to well fit the prior (usually by estimating prior parameters). In the low-intensity case where we generally lack such "useful signals", the final estimation can suffer from a large bias produced by an improperly fitted model. Another previous important contribution to this field is the wavelet-domain hypothesis testing framework first introduced in [START_REF] Kolaczyk | Estimation of intensities of burst-like Poisson processes using Haar wavelets[END_REF][START_REF] Kolaczyk | Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics[END_REF]. Thresholds based on user-specified false detection rate are derived for Haar coefficients both in the constant and the model-based background situations. Universal thresholds are also found available in the above two cases [START_REF] Kolaczyk | Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics[END_REF][START_REF] Charles | Wavelet denoising of Poisson-distributed data and applications[END_REF]. However, this method is only adapted to piecewise constant and burst like intensities as it uses the Haar wavelet, which will yield stair-case-like artifacts in estimating regular intensities. Although Kolaczyk [16] derived the Poisson-corrected version of the Gaussian-based threshold for any wavelet, the asymptotic approximation used in [START_REF] Kolaczyk | Wavelet shrinkage estimation of certain Poisson intensity signals using corrected thresholds[END_REF] may not allow reasonable threshold solution in very low intensity settings. A recent work [START_REF] Zhang | Poisson intensity estimation based on the wavelet domain hypothesis testing[END_REF] has tackled these two drawbacks by proposing a more sensitive bi-orthogonal Haar domain hypothesis testing procedure.

In this paper, we propose a VST for approximately Gaussianizing and stabilizing the variance of a sequence of independent Poisson random variables (RVs) filtered by a lowpass linear filter. This approach is shown to be fast, very well adapted to extremely low-count situations and easily applicable to any dimensional data. We then adapt the VST to a Multi-Scale (MSVST) context to stabilize the detail coefficients of the IUWT of multi-dimensional Poisson count data. We use a hypothesis testing framework in the wavelet domain to threshold the Gaussianized and stabilized coefficients, and then apply the inverse MSVST to get the estimated intensities underlying the Poisson data. Global statistical error rate is also controlled by using a multi-test scheme, i.e. the False Discovery Rate (FDR) procedure.

VST OF FILTERED POISSON DATA

Let (Y n ) be a sequence of RVs observed at the output of a finite impulse response (FIR) filter h ∈ l 2 (Z):

yn = X i h[i]xn-i (1)
where (X n ) are independent Poisson RVs 1 . We additionally assume that the h[i]'s are non-negative samples for i = 1, . . . , N h (typically h is a low-pass filter). Our goal is to stabilize the variance of Y n .

It is known that if h = δ, the Anscombe transform [1] of Y n (hence X n ) acts as if the data arose from a Gaussian distribution with unit variance, under the assumption that the intensity λ n is large. This is why the Anscombe VST performs poorly in low-count settings. But, if the filter h acts as an "averaging" kernel (more generally a low-pass filter), one can reasonably expect that stabilizing Y n would be more beneficial, since the signal-to-noise ratio measured at the output of h is expected to be higher.

A transformation of Y n is sought, such that its variance is constant irrespective of the value of Y n . The general form of the VST A is derived by delta-method argument, giving that:

AY = Z(Y ) = b √ Y + c (2)
In the following, although the general case can be treated easily, for the sake of simplicity we shall assume that ∀n, λ n = λ. For instance, in the a multi-scale wavelet transform context, this amounts to considering that the intensity is constant within the support of the wavelet. The next lemma summarizes the main properties of this transformation:

Lemma 1 (i) From the Taylor series expansion of the RV Z(Y )/b at Y = Y , its mean and variance are given by:

Z/b = √ τ1 √ λ + 4 c τ1 -τ2 8 τ 3/2 1 λ -1/2 + O(λ -3/2 ) (3) Var [Z/b] = τ2 4 τ1 + " 7 τ 2 2 32 τ 3 1 - c τ2 4τ 2 1 - τ3 8 τ 2 1 « λ -1 + " 16c 2 τ2 + 16c τ3 + 5 τ4 64τ 3 1 - 21 c τ 2 2 32τ 4 1 « λ -2 + O(λ -3 ) (4)
where

τ k = i (h[i]) k .
(ii) For the VST to be second order accurate and Z to have asymptotic unit variance, b and c must satisfy:

b = 2 r τ1 τ2 , c = 7τ2 8τ1 - τ3 2τ2 ≥ 0 (5) (iii) For b and c as above, Z -b √ τ 1 λ D → λ→+∞ N (0, 1)
. 1 We can even consider the more general case where X is the sum of mutually independent Gaussian and Poisson variables. But here we prefer to simplify the presentation.

This result tells us that for the chosen value of c, the first order term in the expansion disappears, and the variance is almost constant up to a second order residual term. Thus, for appropriately chosen filter h, this residual is expected to be smaller than for the case of the Anscombe transform applied on the original samples X n . The same reasoning holds for the expansion of Z in (3). Let's consider the case of a 2D B 3 -spline filter 2 (cf. ( 8)), which gives c = 0.0177 and b = 7.3143. With asymptotic unit variance and √ λ for the expectation, the coefficient before the term O(λ -2 ) in ( 4) is 1.72×10 -3 and that before O(λ -1/2 ) in ( 3) is -4.94×10 -4 , while these coefficients associated with the Anscombe VST are respectively -1.56 × 10 -1 and 6.25 × 10 -2 . Therefore, for this kind of filter, the convergence rate towards the asymptotic behavior of Lemma 1 is much faster (about 100 times) for the new VST than for the Anscombe VST. Clearly, the new VST will outperforms the Anscombe VST in low count situations. This is confirmed by the simulations depicted in Fig. 1, where the estimates of Z (resp. Var [Z]) obtained from 500 realizations of Z are plotted as a function of the intensity λ for both Anscombe (dashed-dotted), Fisz (dashed) and our VST (solid). The theoretical bounds from expansions of Lemma 1 limited to the first term (i.e. 1 for the variance and √ λ for the expectation) are also plotted. The saliency of our new VST is obvious. The variance of the RV stabilized using our VST is faster in sticking to the asymptotic bounds. Consequently, the RV transformed using Anscombe can be reasonably considered as stabilized for λ 10, using Fisz for λ 1 and using our VST (after low-pass filtering with h) for λ 0.1.

Fig. 1. Z and Var [Z]

for the Anscombe, Fisz and our (using the B3spline filter) VST. Notice that the RV transformed using Anscombe can be reasonably considered as stabilized for λ 10, using Fisz for λ 1 and using our VST for λ 0.1.

MULTI-SCALE VST AND WAVELET-DOMAIN DENOISING

Although this section focuses on the (not necessarily separable in nD, n ≥ 2) undecimated wavelet transform (UWT), our arguments can be extended similarly to other multi-scale transforms (e.g. the curvelet transform). The UWT W using the filter bank (h, g) of a 1D signal x ∈ ℓ 2 (Z) leads to a set W = {w 1 , . . . , w J , a J } where w j are the wavelet coefficients at scale j and a J are the coefficients at the coarsest resolution.

The passage from one resolution to the next one is obtained using the "à trous" algorithm: aj+1 = h(j) * aj and wj+1 = ḡ(j) * aj

where h[k] = h[-k] and similarly for g. The reconstruction is obtained by:

aj = 1 2 ( h(j) * aj+1 + g(j) * wj+1). (7)
2 Used in the IUWT.

The filter bank (h, g, h, g) needs only verify the exact reconstruction condition. This provides us a high degree of freedom when designing the synthesis prototype filter bank. Because astronomical images contain mostly isotropic sources (stars, galaxies, etc.), astronomers generally prefer using another transform, the Isotropic Undecimated Wavelet Transform (IUWT) [START_REF] Starck | Astronomical Image and Data Analysis[END_REF]. Requirements for a good analysis of such data are: the filters do not need to be orthogonal or bi-orthogonal, but they must be symmetric, and more importantly, h, g, the scaling function φ and the wavelet function ψ must be nearly isotropic. For example, the following filter bank satisfies these requirements:

h 1D = [1, 4, 6, 4, 1]/16 h[k, l] = h 1D [k] • h 1D [l], g = δ -h (8)
From the structure of g, it is easily seen that the wavelet coefficients are obtained just by taking the difference between two resolutions:

wj+1 = aj -aj+1 (9)
In 2D and higher dimensions, at each scale j, we obtain one set {w j } (and not three as in the 2D UWT) which has the same number of samples as the input data. The reconstruction is obtained by a simple co-addition of all wavelet scales and the final smoothed array, namely:

a0 = aJ + J X j=1 wj (10)
That is, the synthesis filters are h = δ and g = δ. Both the analysis and synthesis filter banks implement a frame expansion.

We are now ready to introduce the MSVST-based denoising: a multiscale algorithm for denoising after stabilizing the variance of w j and bringing its distribution closer to normality, when the input signal x is a sequence of independent Poisson RVs. The main steps of this general MSVST denoising algorithm are as follows:

1: Let a 0 = x. For a given filter bank (h, g = δ -h, h = δ, g = δ), 2: for j = 0 to J -1 do 3:
Calculate the approximation coefficients a j+1 using (6).

where A j a j = τ -1/2 1 √ a j + c j is our VST, and c j is the constant in (5) obtained when the scaling function φ(j) plays the role of the low-pass filter considered in the previous section.

5:

Apply the denoising operator D to w j+1 , assuming that they are contaminated by an (almost) zero-mean Gaussian noise, to get the estimates ŵj+1 = D w j+1 . 6: end for 7: Reconstruct â0 :

â0 = A -1 0 AJ aJ + J X j=1 ŵj ! . (12)

Denoising step

This subsection is devoted to some details on the denoising step of the algorithm above. We aim at designing a hypothesis testing-based denoiser in the same vein as in [START_REF] Zhang | Poisson intensity estimation based on the wavelet domain hypothesis testing[END_REF]. To do so, one must access the distributional properties of the stabilized detail coefficients w j , under the null hypothesis that the intensity λ underlying the Poisson process is constant. Thus, Let us consider the RV w j under the null hypothesis. The following proposition gives the asymptotic expansions of w j and Var [w j ]. Proposition 1 We have the following asymptotic expansions:

wj = 0 + O(λ -1/2 ) Var [wj] = 1 τ (j) 1 b 2 j + 1 τ (j+1) 1 b 2 j+1 - ζj+1(h) 2τ (j) 1τ (j+1) 1 + O(λ -1/2 )
where ζ j+1 (h) = k,l H j+1 (k, l)H j (k, l), and:

H j+1 (k, l) = " H j * h(j) (k) h(j) (l) " (k, l), H 0 = δ (13)
where b j (resp. b j+1 ) is defined as in ( 5) associated to φ(j) (resp. φ(j+1) ) and similar for τ . The results above clearly prove that, to a good approximation, the Gaussianized version of the detail coefficients have a zero mean and a variance that depends only on the scale and the chosen filter bank.

Thus, exploiting these results, the hypothesis testing-based denoiser only requires pre-calculating these variances before applying the MSVST, which is computationally simple and fast. Moreover, as the wavelet coefficients are tested simultaneously (multiple testing) and are dependent, we used the FDR strategy under dependency [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] to control the global statistical error.

RESULTS AND DISCUSSION

For illustrative purposes, we have simulated an image with circle-like X-ray sources on a constant background for XMM-Newton telescope. This image can be seen as a model for celestial objects of different size and flux. 3 Each source along any radial branch has the same flux and has a more and more extended support as we go farther from the center. The flux reduces as the branches turn in the clockwise direction. Denoising such an image is highly challenging. The observed image of counts is Fig. 2(a), the restored intensity using our approach in (d) and those restored using some of our competitors Anscombe in (b) and Fisz in (c). For all the methods, the FDR level was fixed at 0.2. Cycle-spinning was used for the Fisz transform in order to not bias the comparison of estimators in favor of our approach.

As revealed by this figure, all estimators perform comparatively well at high intensity levels (image center). However, compared with (b), the relative merits (sensitivity) of the MSVST estimator become increasingly salient as we go farther from the center, and as the branches turn clockwise, i.e. as the intensity becomes low. Even the sources of very low counts were detected by our estimator (see the last branches clockwise in (d) and compare to (b)). Fisz method also proves relatively sensitive. However, it still exhibits a clear staircasing artifact despite the cycle-spinning. Cycle-spinning has the drawback of making the Fisz approach much slower than our method (typically 10 -50 shifts are used in cycle-spinning for each axis of an 512 × 512 image, hence 100 -2500 times slower).

CONCLUSION

In this paper, a multi-scale VST is proposed for intensity estimation of a sequence of independent Poisson RVs. A general wavelet-domain denoiser was also described and has proven very efficient in estimating Poisson noise contaminated data. The algorithm can be easily extended to any dimension and can also be generalized to other multi-scale transforms such as those integrating directionality, e.g. curvelet transform. Another potential application of the MSVST method is the denoising of multi-spectral data, e.g. 2D + energy. Our current research is focusing on these aspects.

=

  Aj aj -Aj+1 aj+1

Fig. 2 .

 2 Fig. 2. (a) Noisy; (b) Anscombe VST; (c) Fisz approach [3]; (d) Our method.

Defined as the integral of the source intensity over its support.


Sylvie Boldo
email: sylvie.boldo@inria.fr

Guillaume Melquiond
email: guillaume.melquiond@inria.fr

Flocq: A Unified Library for Proving Floating-point Algorithms in Coq

Keywords: Floating-point arithmetic, formal proof system, program verification

Several formalizations of floating-point arithmetic have been designed for the Coq system, a generic proof assistant. Their different purposes have favored some specific applications: program verification, high-level properties, automation. Based on our experience using and/or developing these libraries, we have built a new system that is meant to encompass the other ones in a unified framework. It offers a multi-radix and multi-precision formalization for various floating-and fixed-point formats. This fresh setting has been the occasion for reevaluating known properties and generalizing them. This paper presents design decisions and examples of theorems from the Flocq system: a library easy to use, suitable for automation yet high-level and generic.

I. INTRODUCTION

Modern floating-point arithmetic has been widely available for the past thirty years and is by far the most used approach for approximating computations on real numbers. While making it an efficient arithmetic, limits on the precision and exponent range of floating-point numbers are often a cause of unexpected behaviors. Testing of algorithms, software, and hardware, alleviate this issue. Unfortunately, exhaustive testing is usually out of reach and there is no guarantee that all the troublesome inputs have been uncovered.

A higher level of safety can be reached by formal methods in addition to testing. In the broader sense, formal methods cover any approach that provides a mathematical certificate of correctness. The mathematical foundations for such a certificate are given by the IEEE-754 standard [1] which precisely describes the formats and operations. Among formal methods, we find model checking, satisfiability, temporal logic, abstract interpretation, and so on. These approaches are automated and mostly useful on big applications but simple from a numerical point of view.

We are mostly interested in small numerical devices (be they algorithms, programs, microcodes, and so on) that implement state-of-the-art results. Fully automated approaches are useless in this case, so we are left with interactive theorem provers, such as ACL2, Coq, HOL Light, PVS, and so on. They allow the user to guide the reasoning manually, but only experts are able to prove useful properties in a reasonable amount of time. It is therefore important to provide simple yet powerful This work was funded by the F st (ANR-08-BLAN-0246-01) project of the French national research organization (ANR).

formalisms to promote the usage of formal methods and to improve device safety. Bigger systems may not be amenable to a full formal proof but this approach can be applied for critical parts or test cases.

We can distinguish several motivations for the many existing formalizations of floating-point arithmetic. Some of them focus only on formally describing the content of the IEEE-754 and -854 standards for Z [START_REF] Barrett | Formal methods applied to a floating-point system[END_REF], for PVS [START_REF] Carreño | Specification of the IEEE-854 floatingpoint standard in HOL and PVS[END_REF], for Coq [START_REF] Loiseleur | Formalisation en Coq de la norme IEEE-754 sur l'arithmétique à virgule flottante[END_REF]. Others are intended for low-level implementations, e.g. guaranteeing that a circuit or microcode implements a given arithmetic operator for ACL2 [START_REF] Moore | A mechanically checked proof of the correctness of the kernel of the AMD5K86 floating-point division algorithm[END_REF], [START_REF] Russinoff | A case study in formal verification of register-transfer logic with ACL2: The floating point adder of the AMD Athlon processor[END_REF], for Forte [START_REF] Kaivola | Proof engineering in the large: formal verification of Pentium 4 floating-point divider[END_REF], for PVS [START_REF] Jacobi | Formal verification of the VAMP floating point unit[END_REF]. Finally, there are some formalizations dedicated to higher-level properties and algorithms for HOL Light [START_REF] Harrison | A machine-checked theory of floating point arithmetic[END_REF], [START_REF]Formal verification of floating-point trigonometric functions[END_REF], for Coq [START_REF] Daumas | A generic library of floatingpoint numbers and its application to exact computing[END_REF], [START_REF] Boldo | Pitfalls of a full floating-point proof: example on the formal proof of the Veltkamp/Dekker algorithms[END_REF]. To ease the burden of these kinds of proofs, some libraries have focused on automation for Coq [START_REF] Melquiond | Floating-point arithmetic in the Coq system[END_REF], [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF].

We are mainly interested in a formalization of floating-point arithmetic for the Coq proof assistant. As detailed in Section II, there are already several of them with various motivations. As long as the focus is on a particular domain, one can find a suitable library among the existing ones. Unfortunately, the process of formally certifying programs has shown that it was useful to mix several formalisms [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF]. For instance, automation is necessary for dealing with non-overflow proofs and propagation of rounding errors, while a rich formalism is needed for expressing subtle properties such as error-free computations. The usage of several formalisms in a single proof is painful and time-consuming, since it requires one to prove equivalence lemmas between them and to invoke these lemmas repeatedly. That was our motivation: a single library that encompasses several aspects of the existing formalisms.

We have developed the Flocq library available at http://flocq.gforge.inria.fr/ It provides a generic framework: multiple radixes (2, 10, but also odd radixes), multiple precision, several formats for both fixed-point arithmetic and floating-point arithmetic, various rounding modes, and so on (Section III). Special care was taken so that nonstandard arithmetics like flush-to-zero were natively supported. On top of this core library, we have built computable operators to ease automation (Section IV-A). We have also proved several high-level theorems (Section IV-B) that are sometimes improved versions of the theorems found in the existing libraries where they were proved for a single family of formats only.

II. PREVIOUS LIBRARIES IN COQ

The Flocq library is based on our experience using and/or developing floating-point libraries for Coq. This section will describe three of our main sources of inspiration.

A. PFF

PFF is a Coq library initially developed in [START_REF] Daumas | A generic library of floatingpoint numbers and its application to exact computing[END_REF] and various results have been added since then [START_REF] Boldo | Representable correcting terms for possibly underflowing floating point operations[END_REF], [START_REF]A simple test qualifying the accuracy of Horner's rule for polynomials[END_REF], [START_REF] Boldo | Pitfalls of a full floating-point proof: example on the formal proof of the Veltkamp/Dekker algorithms[END_REF]. It has been purely designed to be a high-level formalization of IEEE-754: only floating-point formats with gradual underflow are supported. This is expressed by a formalization where floatingpoint numbers are pairs (n, e) associated with real values n × β e . The requirements for a number to be in the format (e min , β p ) are:

|n| < β p ∧ e min ≤ e.
As there may be several bounded floating-point numbers with the same value, PFF defines a canonical representative when needed. Still, this cohort of equal representatives is sometimes useful as in [START_REF]A simple test qualifying the accuracy of Horner's rule for polynomials[END_REF].

A first advantage of this library is its genericity: the radix is any integer greater than 1 and the precision is any positive integer (it may be 1 in some theorems). It is designed to be as generic as possible in that respect.

As for rounding, the choice was to use an axiomatic approach: there is no function that computes a rounded value, but there is a relation between a real value and a floatingpoint number that says that this floating-point number is a rounding down (up, to nearest, to nearest ties to even. . . ) of the real value. It is practical for proving high-level properties. For example, we can prove theorems for any rounding to nearest, whatever the tie: we simply state that the floating-point number is one of the possibly two roundings to nearest of the real value. A drawback is the lack of automation: to prove that a floating-point number is correctly rounded is very long and tedious.

This formalization is efficient for high-level algorithms containing floating-point technical points. It is convenient for human interactive proofs as shown by the many proofs based on it. The high number of lemmas (about 1400) makes it suitable for a large range of applications.

B. Gappa

Gappa is a tool for automatically proving mathematical properties related to numerical codes [START_REF] Daumas | Certification of bounds on expressions involving rounded operators[END_REF], [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF]. These properties deal with real numbers, but they can tackle floating-and fixed-point arithmetic through the usage of rounding operators, which are functions from R to R. Indeed, the IEEE-754 standard states that "each of the computational operations that return a numeric result specified by this standard shall be performed as if it first produced an intermediate result correct to infinite precision and with unbounded range, and then rounded that intermediate result, if necessary, to fit in the destination's format" [1]. So the result of a floating-point adder, assuming it is not an exceptional value, is just the combination of the real addition with a rounding operator.

When Gappa proves that a property is true, it also generates a formal proof of it. This certificate can be mechanically checked by Coq, and it relies on a support library that contains all the basic facts, e.g. the result of a rounding to binary double precision can be represented with a mantissa of at most 53 bits. Gappa has been designed for radix-2 arithmetic; therefore, that is the only radix the support library has theorems for.

So that generated certificates can be checked without any human intervention, they rely on computations performed inside the proof checker. In particular, given a number m • 2 e , the support library is able to correctly round it to a number m • 2 e that fits in a given format. The library is not restricted to some specific formats; it can handle any floating-or fixedpoint format that can be described by an integer relation between exponents. This mechanism has been reused for Flocq and is described in more details in Section III-D.

To summarize, Gappa's support library is a complete formalization of binary arithmetics. It also provides some computable operators for performing rounding. But it is designed for a usage by proofs generated automatically, so its theorems have awkward statements that are not really suitable for interactive proofs.

C. Coq.Interval

The Coq.Interval library [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF] provides a tool based on interval arithmetic and Taylor series for automatically proving inequalities on real-valued expressions inside Coq. The methods are well-known, if not outdated, but their primary characteristic is that they are performed inside the logic formalism of the proof checker. The key concept is that the whole formalization has to be computable, in order for the interval operators to return actual results. They rely on a computable floating-point arithmetic that is formalized in Coq [START_REF] Melquiond | Floating-point arithmetic in the Coq system[END_REF]. It supports a single family of formats: multiprecision with unlimited exponent range. From an interface point of view, the library is therefore similar to MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF], but it supports any radix.

Unlike Gappa, there is no external oracle for precomputing arithmetic results this time. So Coq.Interval not only provides computable rounding operators but also any operator the interval methods may need: division, square root, some elementary functions, and so on. Some of these algorithms have been reimplemented in Flocq and are described in Section IV-A.

Since this floating-point kernel is primarily intended for interval arithmetic, rounding modes other than directed toward infinities have been neglected in the formalization. Moreover, the library does not provide many facts about floating-point arithmetic, except its own proofs of correctness.

III. CORE LIBRARY

In this section, we will describe the core of the library, namely the main definitions and choices of representation.

A. Rounding predicates

Here, we consider the floating-point numbers as any subset F of the real numbers. This means there is no requirement that F be discrete or that 0 ∈ F. There will be two definitions of the rounding modes. The first one, described here, is that of Section II-A: rounding modes are relations between real numbers and their rounded values. The second definition, based on computable functions, will be described in Section III-E but it does not apply to the most exotic formats.

A rounding predicate Q has the Coq type R → R → Prop, meaning it is a relation between two real numbers (as ≤ or =). For it to be a rounding predicate, it must fulfill two requirements. First, it must be total (any real can be rounded): ∀x ∈ R, ∃f ∈ R, Q(x, f ). Second, it must be monotone: ∀x, y, f, g ∈ R, Q(x, f ) ⇒ Q(y, g) ⇒ x ≤ y ⇒ f ≤ g (that is to say nondecreasing). These two properties are enough for a rounding predicate to have reasonable properties. For example, those properties imply the uniqueness of rounding: two reals that are roundings of the same value are equal.

A more interesting property is the fact that, from any rounding predicate Q, we can construct a function r such that for all real x, Q(x, r(x)).

We can then define the common rounding modes by their mathematical properties as in Section II-A for a given F:

• Rounding toward -∞ is denoted by (x, f ) and defined by f ∈ F ∧ f ≤ x ∧ (∀g ∈ F, g ≤ x ⇒ g ≤ f ). • Rounding toward +∞ is denoted by (x, f ) and defined by f ∈ F ∧ x ≤ f ∧ (∀g ∈ F, x ≤ g ⇒ f ≤ g). • Rounding toward 0 is defined by (0 ≤ x ∧ (x, f )) ∨ (x ≤ 0 ∧ (x, f )).
It is easy to prove that these directed modes are monotone and that a rounded value is unique (if (x, f 1 ) and (x,

f 2 ) then f 1 = f 2 )
. But what about the existence of these values?

For instance, if 0 is not in F, strange things happen: consider for example F = {n×2 emin with n ∈ Z * } in radix-2. Then 0 ∈ F and if we round down the positive value 2 emin-1 , we get the negative value -2 emin , and -2 emin-1 is rounded up to 2 emin . So such a rounding toward zero would not be monotone.

The set F must also be such that we are able to round down: ∀x ∈ R, ∃f ∈ R, (x, f ). This is needed as roundings do not exist for any F: consider for example the set of rational numbers. Then an irrational number x cannot be rounded: there is always another rational that is smaller than x and nearer to x.

At last, we assume the symmetry of F. This is not a necessary condition, but it is reasonable and greatly helps the following proofs. The idea is that, if we know how to round down, a symmetric set allows us to round up for free:

(x, f ) ⇔ (-x, -f ).
Theorem 1 (satisfies any imp {DN,UP,ZR}) For any set F such that 0 ∈ F, F is symmetrical, and rounding to -∞ exists, then rounding toward -∞, toward +∞, and toward zero are rounding predicates (they are total and monotone).

B. Rounding to nearest

The most used (and default mode) is rounding to nearest, ties to even (when in the middle, choose the one with the even mantissa). But the revision of the IEEE-754 standard [1] added the rounding to nearest, ties to away (when in the middle, choose the one with the biggest absolute value). With these two roundings having much in common, it is natural to define a generic rounding to nearest, whatever the tie-breaking rule. So we first define rounding to nearest without ties by f ∈ F∧

(∀g ∈ F, |f -x| ≤ |g -x|) and denote it by •(x, f ).
This rounding to nearest does not have a tie-breaking rule: it is not equivalent to a function as two different real numbers may be correct roundings to nearest of the same real (when it is the middle of two floating-point values). For ties away from zero, the definition is easy:

•(x, f )∧(∀g, •(x, g) ⇒ |g| ≤ |f |).
For a more generic tie-breaking rule, we define a rounding to nearest that depends on a property P . This generic rounding to nearest with ties is denoted by • P (x, f ) and defined by

•(x, f ) ∧ (P (x, f ) ∨ ∀g, •(x, g) ⇒ f = g).
This means that either f has the property P , or it is the only rounding to nearest. You may think of P as being "f is even" or |f | ≥ |x|. This generic rounding is unique and monotone under the following assumption on P : if (x, d) and (x, u) and • P (x, d) and • P (x, u) and P (x, d) and P (x, u), then we have d = u. We also proved that the symmetry of P implies the symmetry of • P and that rounding to nearest, ties away from zero is exactly this generic rounding to nearest with

P (x, f ) def = |f | ≥ |x|.
The most difficult point is of course the existence of rounding to nearest. As far as rounding to nearest without ties is concerned, the existence is guaranteed as soon as the three properties of the preceding subsection are fulfilled (but it is not a rounding mode as it is not monotone). For the generic rounding to nearest, we need this property on P : if x ∈ F and (x, d) and (x, u), then we either have P (x, d) or P (x, u). If this last property (and the three preceding ones) are fulfilled, then • P is a rounding mode. We easily deduce that the rounding to nearest, ties away from zero is a rounding mode.

For rounding to nearest ties to even, we cannot define it here: the reason is that we cannot decide if a value is even or odd. Let us consider x = 2 e ∈ F: is it odd as x = 1 × 2 e or is it even as x = 2 × 2 e-1 ? We need both a unique canonical representation of a floating-point number and a definition of ties to even. The representation will be defined in Section III-D and the rounding will be defined in Section III-E.

C. Floating-point numbers

Now, let us define the most useful formats. They will be defined differently in the next section, but we begin with the most intuitive definitions. We consider floating-point numbers as a mantissa and an exponent. The mantissa is shifted in order to get an integer rather than a fixed-point value. A floatingpoint number f is then a pair of integers (n f , e f ) and its value is F2R(f ) = n f • β e f , where β, the radix, is an integer greater than one.

From that, we can define a format as a set of real numbers corresponding to all the expected floating-point numbers. For example, the set of fixed-point numbers with exponent e min is the set of real numbers such that there exists a floatingpoint number with exponent e min having this value. Note that a format is therefore a subset of R such that there exists a floating-point number equal to the real that has such and such property.

Here are the common formats: fixed-point (FIX), floatingpoint with unbounded exponents (FLX), normalized floatingpoint with unbounded exponents (FLXN), floating-point with gradual underflow (FLT), and floating-point with flush-to-zero (FTZ). The array below formally defines these formats.

Format is defined by ∃f, F2R(f ) = x ∧ . . . FIX emin (x) e f = e min FLX p (x) |n f | < β p FLXN p (x) x = 0 ⇒ β p-1 ≤ |n f | < β p FLT p,emin (x) e min ≤ e f ∧ |n f | < β p FTZ p,emin (x) x = 0 ⇒ e min ≤ e f ∧ β p-1 ≤ |n f | < β p
To ease further proofs, we proved that the FLX and FLXN formats are equivalent.

These are high-level definitions, and the FLT format is exactly the format used in PFF (see Section II-A). Nevertheless, these are not the definitions we will mainly use: we will base the format on the definition of a ϕ function and prove the equivalence between both definitions (see Section III-D).

D. Generic format

Representing formats as predicates and rounding modes as relations on real numbers makes it simple to manipulate them when rounded values are known beforehand. But they are a hindrance when a proof actually needs to compute a rounded value. Therefore we have chosen a more computational approach for another representation of formats and rounding modes.

In order to get closed formulas for rounded values, we had to put constraints on the number sets. All the formats of this family (called generic format) satisfy the following two main properties: they contain only floating-point numbers m • β e and all the representable numbers in a given slice are equally distributed.

The slice of a real number x is given by its discrete βlogarithm e = slice(x):

β e-1 ≤ |x| < β e .
F ϕ is entirely described by a function ϕ : Z → Z that transforms numbers' discrete logarithms into canonical exponents for this format. In other words, a number x is in F ϕ if and only if it can represented as a floating-point number m • β ϕ(slice(x)) . Figure 1 shows the graph of the canonical exponents for the three usual families of floating-point formats.

More precisely, the canonical exponent of a real number x is cexp(x) = ϕ(slice(x)). Its scaled mantissa is smant(x) =

x • β -cexp(x) . And the format F ϕ is defined as the set of real numbers x such that

x = Z(smant(x)) • β cexp(x)
with Z the integer part (rounded toward zero).

The above definition is equivalent to saying that smant(x) is exactly an integer; but stating it as a rewriting rule makes it a bit simpler to use in Coq. From there, one can deduce that generic formats contain zero and are symmetric. A slice e contains representable numbers, e.g. β e-1 , when ϕ(e) < e.

In order to show that F ϕ is a suitable format as described in Section III-A, there is a property left to satisfy: any real number should have a rounded-down value. This value can be chosen as (x) = smant(x) • β cexp(x) , but we have to prove that this function is both increasing and onto F ϕ . Note that the function is the identity for values of F ϕ , by definition of F ϕ . Similarly, the rounded-up value would ideally be (x) = smant(x) • β cexp (x) . What are the constraints on function ϕ so that these functions are proper rounding?

First, let us consider the slice [β e-1 , β e ). Let us assume that ϕ(e) < e, so that the slice contains some representable numbers. Ideally, a real number x = β e -ε (with ε positive but negligible) should be rounded up to (x) = β e , but this power might not be representable since it lies in a different slice (e + 1). In order to ensure that it is representable, we put the following constraint on ϕ: ϕ(e) < e ⇒ ϕ(e + 1) ≤ e.

Second, let us consider a real number x inside a slice e such that e ≤ ϕ(e). We have (x) = 0 and (x) = β ϕ(e) . Zero is representable, but β ϕ(e) might not be. Moreover, all the other real numbers inside the open interval (0, β ϕ(e) ) have to round down to 0 and up to β ϕ(e) too; otherwise the two functions would not be increasing. This requires the following constraint on ϕ: e ≤ ϕ(e) ⇒ ϕ(ϕ(e) + 1) ≤ ϕ(e), ∀e , e ≤ ϕ(e) ⇒ ϕ(e ) = ϕ(e).

Above constraints are necessary; but we have also proved that they are sufficient for and to be rounding functions on F ϕ . If ϕ satisfies these constraints, it is said to be valid (predicate valid exp). This is true for all the usual formats and we will give the corresponding ϕ functions below. The formats described by these functions are proved to be equivalent to the formats of Section III-C.

A fixed-point format containing all the multiples of β emin is described by FIX exp(e) = e min .

An unbounded floating-point format with a precision of p digits is described by The ulp (unit in the last place) of a real number x is defined as β cexp(x) . This function is partial: it is not defined for 0 (neither was cexp). But it is defined for any other real, be it in F ϕ or not. An immediate property of ulp is Theorem 2 (ulp DN UP) Assuming that ϕ is valid, ∀x ∈ F ϕ , (x) = (x) + ulp(x).

A generic format is said to not have the flush-to-zero property if ulp(x) is representable in F ϕ for any real number x. The corresponding property on ϕ is predicate not FTZ prop: ∀e ∈ Z, ϕ(ϕ(e + 1)) ≤ ϕ(e).

E. Rounding operators for generic formats

The expressions of the and rounding operators above can be generalized by using a function Zrnd : R → Z:

round ϕ (x) = Zrnd(smant(x)) • β cexp(x) .
For round ϕ to be a rounding mode, Zrnd has to be the identity for integer inputs and to be increasing on R. That way, rounding modes are given by closed formulas instead of relations. Due to the monotonicity requirement, this is more than being a faithful rounding as defined by Priest [START_REF] Priest | Algorithms for arbitrary precision floating point arithmetic[END_REF].

The expression above is not the most generic, since Zrnd(x) has only access to the scaled mantissa and not the original real number. As a consequence, round ϕ will behave in a similar way for all the slices that have the same precision. An additional argument could have been used for Zrnd, e.g. the canonical exponent, but we did not find a reasonable setting that was worth the noise.

As shown previously, choosing Zrnd(m) = m (resp. m ) gives an operator that rounds downward (resp. upward). We have proved that they are equivalent to the rounding relations of Section III-A.

Rounding to nearest (with tie breaking to even) is not much harder to define: the scaled mantissa m just has to be rounded to the nearest integer, tie breaking to the even one. While the definition is simple, one still has to prove that the resulting operator satisfies the expected property on parity. Indeed, even if the integer is even, the scaled mantissa of the rounded value might be odd if it lies in a different slice. To avoid this parity discrepancy, either the radix should be odd, or the following properties should hold in addition to the previous constraints on ϕ:

ϕ(e) < e ⇒ ϕ(e + 1) < e, e ≤ ϕ(e) ⇒ ϕ(ϕ(e) + 1) = ϕ(e).

We have proved that, under these requirements, the computable rounding to nearest is equivalent to a rounding to nearest that is the one described in Section III-B with a P that is "the canonical floating-point number (with exponent cexp(x)) has an even mantissa".

These properties on ϕ are satisfied by the FIX format. For FLX and FLT formats of precision p, these properties degenerate to simply 1 < p. For p = 1, we have neither the existence, nor the uniqueness of rounding to nearest, ties to even. The case of the FTZ format is special, since breaking tie to even is no longer deterministic: both 0 and its successor β p-1 • β emin have even mantissas when 1 < p and β is even.

Note that rounding operators do not imply that we cannot use a generic rounding to nearest. We defined a Zrnd to nearest depending on a choice function that gives either the rounding to nearest when unique or one of the two possible values depending on the result of the choice function when in the middle. By proving some theorems whatever the choice function, we retrieve the behavior of the generic rounding to nearest without ties defined in PFF and in Section III-B.

Another interesting property of format FTZ is the way the rounding operators are defined on the subnormal range. There are two possible implementations, either small real values are correctly rounded toward the smallest normal numbers (sometimes called abrupt underflow), or they are simply flushed to zero, irrespective of the rounding direction. In the first case, the previous definitions of the rounding operators are still valid. In the second case, a variant must be used:

Zrnd FTZ (x) = Zrnd(x) if |x| ≥ 1, 0 otherwise.

IV. OTHER PARTS A. Computable functions

There are some situations where an axiomatic approach to floating-point arithmetic is no longer sufficient and one would like to actually compute with a proof assistant. That may be useful for generating test values, or for brute-forcing the few inputs that elude the correctness proof of an algorithm [START_REF] Harrison | Formal verification of square root algorithms[END_REF], or for automating proofs of theorems, be they about floatingpoint arithmetic [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF] or real arithmetic [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF].

1) Computable rounding operators: While the rounding formulas of Section III-E are closed, they be used to effectively compute with real numbers. A subset on which

• -like functions are decidable has to be chosen. This could be the subset of rational numbers; it would be suitable for dealing with addition, multiplication, and division of floatingpoint numbers, but difficulties would arise for square root already. We have chosen a different approach: we have decided to over-approximate real numbers by a floating-point number (not necessarily in F ϕ ) and a location relative to this number. Let us consider a real number x and a floating-point number d = m•β e such that m•β e ≤ x < (m+1)•β e = u. This is the property inbetween float which relates x, m, e, and a location . If x is equal to d, then tells so: = loc Exact. Otherwise, tells whether x is smaller, equal, or bigger than (d + u)/2. So there are four different locations of x with respect to d. They are similar to the usage of the round and sticky bits that can be found in hardware implementations.

This approach restricts the rounding operators that can be computed by our formalism. They have to round all the real numbers strictly between d and (d + u)/2 toward the same floating-point number of F ϕ , and similarly for (d+u)/2 and u, assuming that d and u are consecutive floating-point numbers in F ϕ . Fortunately, this property is sufficient to express all the standard rounding modes, as shown by Theorems inbetween float DN, UP, NE, and so on.

Given a lower bound d = m • β e and a location of the real number x, one can compute a new lower bound d = m/β • β e+1 and a location of x with respect to d . This location is computed from β, the remainder of the Euclidean division of m by β, and only. The process can be iterated until one obtains a lower bound in F ϕ . The resulting location then tells how to round the real number.

More precisely, cexp(d) = ϕ(slice(m)+e) can be computed by counting the digits of m and applying ϕ. If x ≥ 0 and cexp(d) ≤ e, then cexp(x) = cexp(d). A suitable triple (m , cexp(x), ) (Theorem truncate correct) can then be computed. At this point, Theorem round any correct explains how to perform the final rounding. Some simpler theorem instances that combine both theorems with a given rounding mode are provided for ease-of-use: Note that the computed floating-point number does not necessarily have a canonical exponent.

2) Computable arithmetic operators: Being able to compute a correctly-rounded floating-point number from a triple is only part of the work. One should also be able to compute triples for basic arithmetic operations. For addition and multiplication of floating-point numbers, the result of the real operation is representable as a floating-point number, so building a suitable triple is straightforward.

The situation is slightly more complicated for division and square root, since their results are often not representable. As a consequence, the operators Fdiv core and Fsqrt core also take a precision as input. The triple they compute is then guaranteed to locate the real result and to have at least the requested number of digits. Both operators are similar: they first shift the mantissas of the floating-point inputs and then perform an integer division or square root. The integer remainder is used to decide the location of the real result with respect to the computed floating-point number. The implementation of these operators was taken from the floatingpoint kernel of Coq.Interval [START_REF] Melquiond | Floating-point arithmetic in the Coq system[END_REF].

Because the operators compute only triples and not rounded results, they are independent of formats and rounding modes (as can be seen by the position of Fcalc div and Fcalc sqrt on Figure 2). It it up to the user to combine them with rounding operators to get a full rounded operator. This combination amounts to choosing the number of digits of the triple and proving that it is an upper bound on the size of the rounded result. For instance, always asking for p digits when the format is FLT with precision p is sufficient, whether the result is normal or subnormal. The core of the correctness proof of such an operator is about 15 lines, which makes it tractable from a Coq point of view. The library provides Theorem Fsqrt FLT ne correct as an example of instantiating the square root operator for a specific format and a specific rounding mode.

Note that these arithmetic operators were designed with correctness and computability in mind, not performance. For instance, adding to x a negligible number y will cause the adder to shift x until it has the same exponent as y, then the result will be truncated until x (or some value close to it) is obtained again. The efficient way of performing the addition would have been to truncate y until it has the same exponent as x instead. For other operators, performance issues will arise when the inputs have a precision much bigger than the output. This is nothing new; efficient algorithms have already been designed for the multi-precision library MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF] but their correctness proofs are on a different level than the proofs of Flocq's operators.

B. High-level lemmas

To assess the usefulness and practicality of our library, we have chosen to prove well-known generic-radix theorems that are either helpful, or expectedly difficult, meaning that they are proved in PFF but they belong to where the formalizations differ the most.

1) Exact subtraction: This well-known fact [START_REF] Sterbenz | Floating-Point Computation[END_REF] states that if x and y are floating-point numbers near enough one to another (y/2 ≤ x ≤ 2y), then the subtraction x-y is computed without rounding.

Theorem 4 (sterbenz) Assuming that ϕ is valid satisfies ∀e x , e y ∈ Z, e x ≤ e y ⇒ ϕ(e x ) ≤ ϕ(e y ), Then, for all x and y in F ϕ such that y 2 ≤ x ≤ 2y, we have that x -y ∈ F ϕ .

As our format is very generic, there is a requirement on the ϕ function, that is that ϕ is monotone. This is reasonable as both FIX, FLX, and FLT fulfill it. This is also necessary: let us consider in radix 2 a function ϕ such that ϕ(0) = -2 and ϕ(-1) = -1. Let x = 0.75 and y = 0.5. The slice of both x and y is 0 as 2 -1 ≤ x, y < 2 0 . So x and y are in the generic format with exponent ϕ(0) = -2. Moreover, y/2 ≤ x ≤ 2y. Then x -y = 0.25 has a slice of -1 and therefore needs to be represented with an exponent ϕ(-1) = -1 but 0.25 cannot be represented as m × 2 -1 .

2) Relative error: The idea is that there exists a small ε such that round(x) = x • (1 + ε). Then rounded values are replaced with such formulas and the various ε are added, multiplied, and so on. See [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] for examples of use.

We proved this property also with a requirement on the ϕ function:

Theorem 5 (generic relative error ex) Let us assume that ϕ is valid and that there exists p and e min such that

∀k ∈ Z, e min < k ⇒ p ≤ k -ϕ(k).
Then, for any rounding operator round ϕ and for any real x such that β emin ≤ |x|, there exists ε such that |ε| < β 1-p and round ϕ (x) = x • (1 + ε).

We also proved the special case of rounding to nearest where |ε| ≤ β 1-p 2 . Note that the requirement on the ϕ function intuitively corresponds to the non-underflow cases of the FLT format. This cannot be applied to fixed-point formats for example.

From that theorem, we deduced its application to FLX and FLT: in FLX, there is no requirement on x for the theorem to hold with p being the precision. In FLT, we pose e min and p as expected and we have the existence of ε provided β emin+p-1 ≤ |x|.

3) Error of the addition: It is a well-known fact that, using rounding to nearest, the error of a floating-point addition is a floating-point number. It can even be computed using floatingpoint operations [START_REF] Knuth | The Art of Computer Programming[END_REF]. This theorem and the following ones about multiplication, division, and square root are the bases of error-free transformations [START_REF] Ogita | Accurate sum and dot product[END_REF].

As before, we need a hypothesis on the ϕ function (the same as for Sterbenz exact subtraction): Theorem 6 (plus error) Assuming that ϕ is valid and monotone, for any rounding to nearest round N ϕ ,

∀x, y ∈ F ϕ , round N ϕ (x + y) -(x + y) ∈ F ϕ .
The genericity of the theorem with respect to formats can be seen on Figure 2: the Fprop plus error file does not depend on the specialized formats defined in Fcore FIX, FLX, and so on.

4) Error of the multiplication: It is also well-known that the error of a floating-point multiplication is a floating-point number, whatever the rounding, but provided no underflow occur [START_REF] Knuth | The Art of Computer Programming[END_REF].

We then first proved this using the FLX format:

Theorem 7 (mult error FLX) For any rounding operator round ϕ X on a FLX format

F ϕ X , ∀x, y ∈ F ϕ X , round ϕ X (x × y) -x × y ∈ F ϕ X .
To get a result in the FLT format similar to [START_REF] Boldo | Representable correcting terms for possibly underflowing floating point operations[END_REF], we first proved the preceding theorem while giving precisely its exponent. With this fact and results about translation between FLX and FLT, we easily proved the corresponding theorem in FLT with additional hypothesis for the gradual underflow: Theorem 8 (mult error FLT) For any rounding operator round ϕ T on a FLT format F ϕ T with a minimal exponent e min , ∀x, y ∈ F ϕ T , β emin+2p-1 ≤ |x × y| ⇒ round ϕ T (x × y) -x × y ∈ F ϕ T . 5) When the rounding is zero: Another useful fact is when x + y rounds to zero, then it usually means that x + y is mathematically equal to zero. This is correct in FLT and in FLX. We proved it for a generic format: Theorem 9 (round plus eq zero) Assuming that ϕ is valid and satisfies the predicate not FTZ prop: ∀e ∈ Z, ϕ(ϕ(e + 1)) ≤ ϕ(e), for any rounding operator round ϕ , ∀x, y ∈ F ϕ , round ϕ (x + y) = 0 ⇒ x + y = 0.

The assumption on ϕ seems strange. It corresponds to the fact that the format must not be FTZ. Indeed, this theorem does not hold in FTZ and this requirement (that is fulfilled in FLX, FLT, and FIX) is sufficient to avoid the bad cases.

6) Remainder of the division and square root: It is wellknown that the of the division, that is to say the real value x -•(x/y) × y, and the remainder of the square root, that is to say x -•( √ x) × •( √ x), usually are floating-point numbers. It is not the case when an underflow occurs [START_REF]A simple test qualifying the accuracy of Horner's rule for polynomials[END_REF]. This is the reason why we proved it in FLX format, meaning we assume there is no underflow.

Theorem 10 (div error FLX) For any rounding operator round ϕ X on a FLX format F ϕ X , ∀x, y ∈ F ϕ X , x -round ϕ X (x/y) × y ∈ F ϕ X .

The square root requires rounding to nearest.

Theorem 11 (sqrt error FLX) For any rounding to nearest round N ϕ X on a FLX format

F ϕ X , ∀x ∈ F ϕ X , x -round N ϕ X ( √ x) 2 ∈ F ϕ X .
In contrast to intuition, the two preceding theorems do not have to assume anything on x and y (hence are easier to apply), due to some peculiarities of real arithmetic in Coq. For division, y does not have to be nonzero as x-round ϕ X (x/y)×y is provably equal to x when y = 0. The square root of a negative number is defined as zero in Coq, so x -(round N ϕ X ( √ x)) 2 is provably equal to x when x < 0.

7) Predecessor and successor: Floating-point numbers are a discrete set. This is still true for our generic format with no additional requirement on the ϕ function. Among the useful properties of the fact that the set of floating-point numbers is discrete is the existence of a predecessor and a successor. For x > 0 in the format, its successor is x + ulp(x) and its predecessor is

pred(x) = x -β ϕ(slice(x)-1) if x = β slice(x)-1 , x -ulp(x)
otherwise.

These values are proved to be in the format and are such that (x + ε) = x + ulp(x) and (x -ε) = pred(x) (if pred(x) > 0) for a small enough ε.

Theorem 12 (pred ulp) Assuming ϕ is valid, for any x ∈ F ϕ such that 0 < pred(x), we have pred(x) + ulp(pred(x)) = x.

Theorem 13 (le pred lt, succ lt le) Assuming ϕ valid, for any x and y in F ϕ such that 0 < x < y, we have x ≤ pred(y) and x + ulp(x) ≤ y.

V. CONCLUSION

In the formal proof community, people are aware that designing libraries is a challenge, as they should be both usable and correct. For example, [START_REF] Barrett | Formal methods applied to a floating-point system[END_REF] had subtle discrepancies with the IEEE-754 standard, as pointed out by Pr Kahan at the July 2002 meeting of the IEEE 754R group. 1 The correctness of Flocq comes from its mathematical definition of rounding predicates: it is hard to get them wrong when rounding down is defined as the biggest floating-point number smaller than a given real number. As for its usability, it was exercised by the many high-level theorems we were able to prove in a convenient way.

This paper presents only a small overview of the library, as it contains more than 450 theorems. The library is available under an open-source license at http://flocq.gforge. inria.fr. Its definitions and theorems can also be browsed online. The library has 23 files in three directories: the core of the library, its computable definitions and properties, and the high-level properties. Figure 2 shows how the files from the three directories depend on each other.

As a consequence of the generalization of the theorems to other formats (including possibly degenerate ones), this work has also been the occasion to uncover some pathological cases. In particular, our decision of natively handling flushto-zero formats has led us to consider unusual requirements for common theorems. This extreme generalization puts Flocq apart from the existing libraries (be they for Coq or other systems) which tend to focus on one specific arithmetic.

Since Flocq was motivated by our experience with previous libraries (in particular PFF and Gappa), we will conclude by comparing them to Flocq.

A. PFF

The main characteristics of PFF are its practicality for the FLT format, its size, and its lack of automation. As the roundings are axiomatic, it is very difficult to prove that a given floating-point number is the rounding of a given real. Flocq solves this problem by easing the usage of automated tools (see Section V-B).

Compared to the PFF library, Flocq is more generic as it easily handles both FLT, FLX, FIX, and even more exotic formats such as formats without subnormal numbers. As PFF was dedicated to FLT, it has a different basic formalization from Flocq, therefore theorems cannot be automatically translated from one to another. This is unpleasant as it means re-doing all the proofs.

The worst expected case for translation is when the proofs of PFF heavily rely on the representation of floating-point numbers, such as the proofs of the facts that the error of the addition, multiplication, the remainder of the division and of the square root are representable floating-point-numbers (under various assumptions, see Section IV-B). This is the reason we chose them as experiments for the translations of theorems. Those pure FLT proofs have been translated and Fig. 2. Dependency graph of the library. Theorems from one file depend on the theorems from the pointed files. The shapes of the node classify the origin of the files: rectangles for the core library (Section III), octagons for the computable operators (Section IV-A), and ellipses for the high-level properties (Section IV-B).

generalized to more formats. The proofs had to be rethought as the basic formalizations were different. Even if the proof ideas are mainly the same, the details and the hypotheses are different.

Nevertheless, our experiments have proved it is possible. We are currently proving properties about the quasi-correctness of •(•(a × x) + y), when |y| is much greater than |a × x| [START_REF]A simple test qualifying the accuracy of Horner's rule for polynomials[END_REF].

B. Gappa

Two of the motivations for the Flocq library were: avoiding proof duplication between formalisms for Coq and allowing automation within a high-level formalism. Gappa's support library has therefore been rewritten to depend on Flocq.

The first benefit of this change is a reduction of the size of the library by a third, while supporting more theorems. For instance, the usual semantic of rounding modes is now available for Gappa's operators, e.g. (x) is the biggest representable number that is less or equal to x. Such properties were already true before, but they had to be proved outside Gappa's formalism [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF]. More generally, Gappa's rounding operators are now simply Flocq's ones and hence benefit from all their theorems.

The second benefit of Gappa no longer defining its own operators is that the gappa tactic of Coq can be used to automate the proof of goals written with Flocq's formalism. Consider the example of Figure 3. Line 1 defines the standard double-precision floating-point format: radix-2 FLT format with a precision of 53 bits and a minimal exponent equal to -1074. Line 3 is just a notation for rounding to zero in order to shorten the proof. The goal to prove ranges from line 7 to 11; it states that if the real numbers a ∈ [52/16; 53/16] and b ∈ [22/16; 30/16] are representable in the format, their difference is representable too. Note that this is not an instance of exact subtraction as stated in Section IV-B1, since a/b is possibly bigger than 2. The proof of this property is lines 13-18. First, line 14 states that proving that a -b is representable is the same as proving that it is equal to its rounded value (see Section III-D). For the same reason, lines 16 and 17 replace a and b by their rounding to zero. Finally, the Gappa tool is called to automatically complete the proof. The completed proof is entirely checked by Coq.

C. Perspectives

Certifying the proof obligations associated to the correctness of numerical programs requires mixing several floatingpoint formalisms [START_REF] Boldo | Combining Coq and Gappa for certifying floating-point programs[END_REF]. The tediousness of this approach was an incentive for developing the Flocq library. These proof obligations are generated by an automated tool: the Frama-C/Jessie/Why toolchain. It takes an annotated C program as input, performs weakest precondition computations on it, and generates theorem statements that, once proved, guarantee that the original program fulfills its specification. Why was using the formalism of PFF for generating the floating-point obligations [START_REF] Boldo | Formal verification of floating-point programs[END_REF]. It now uses Flocq.

Another objective is to fill the gap between a high-level formalization of floating-point arithmetic and the IEEE-754 description as bit vectors. The goal is to have the tools for a full certification of an algorithm from its C source code to its final assembly code. One of the building blocks is CompCert, a certified compiler written in Coq [START_REF] Leroy | A formally verified compiler back-end[END_REF] and floating-point support is being added to compiler. Therefore, related code constructs and optimizations have to be formally proved too. For instance, the conversion from integers to floating-point values may involve a mix of bit-level operations (writing a magic constant in the most significant word of a floating-point register) and floating-point operations. Proving the correctness of such a code sequence requires an extended formalization.

We have started such a development to handle exceptional values (signed zeros, infinities, NaN) and behaviors.

FLXFig. 1 .

 1 Fig.1. Values of ϕ for formats FLX, FLT, and FTZ, with precision p. These functions are the same for normal numbers (e ≥ e min + p), but they diverge for subnormal numbers. In particular, the ϕ function for FTZ is discontinuous.

Theorem 3 (

 3 round trunc UP correct) Assuming ϕ is valid, for any positive real number x over-approximated by a triple (m, e, ), if either e ≤ ϕ(digits(m) + e) or = loc Exact, then (m , e , ) = truncate ϕ (m, e, ) ⇒ round UP ϕ (x) = m • β e if = loc Exact, (m + 1) • β e otherwise.

1 4 (

 4 Definition format := 2 generic_format radix2 (FLT_exp (-1074) 53). 3 Notation rnd := round radix2 (FLT_exp (-1074) 53) rndZR).

14 change

 14 (a -b = rnd (a -b)).

15 revert

 15 Ia Ib.

16

  replace a with (rnd a).

17

 17

  replace b with (rnd b).

Fig. 3 .

 3 Fig.3. Example of a proof script using automation.

  1 http://grouper.ieee.org/groups/754/meeting-minutes/02-07-18.html

				Fprop_div_sqrt_error	Fprop_plus_error	Fprop_Sterbenz
	Fcalc_sqrt	Fcalc_div	Fcalc_round	Fprop_relative	Fprop_mult_error
		Fcalc_digits	Fcalc_bracket	Fcore	Fcore_FTZ	Fcalc_ops
				Fcore_FLT	Fcore_FLX	Fcore_FIX
						Fcore_ulp
						Fcore_generic_fmt
						Fcore_float_prop	Fcore_rnd
						Fcore_defs	Fcore_Raux