Steven Gazal

Mourad Sahbatou

Marie-Claude Babron

Emmanuelle Génin

Anne- Louise Leutenegger

FSuite: exploiting inbreeding in dense SNP chip and exome data

FSuite is a user friendly pipeline developed for exploiting inbreeding information derived from human genomic data. It can make use of SNP chip or exome data. Compared to other software, the advantage of FSuite is to provide a complete suite of scripts to describe and use the inbreeding information. It includes a module to detect inbred individuals and estimate their inbreeding coefficient, a module to describe the proportion of different mating types in the population and the individual probability to be offspring of different mating types that can be useful for population genetic studies. It also allows the identification of shared regions of homozygosity between affected individuals (homozygosity mapping) that can be used to identify rare recessive mutations involved in monogenic or multifactorial diseases.

INTRODUCTION

Inbreeding is a central concept in genetics. In population studies, the inbreeding coefficient f of individuals is evaluated to characterize mating habits. In rare disease studies, homozygosity mapping [START_REF] Lander | Homozygosity mapping: a way to map human recessive traits with the DNA of inbred children[END_REF] has been widely and successfully used to localize variants with strong recessive effect by searching for regions homozygous by descent (HBD) on pedigrees with inbred cases.

With the availability of dense genome-wide genetic data, it is now possible to study inbreeding for individuals without genealogical information. Several software packages are available to estimate inbreeding coefficients using genetic data. Some of these provide singlepoint estimates of the inbreeding coefficient such as PLINK [START_REF] Purcell | PLINK: a tool set for whole-genome association and populationbased linkage analyses[END_REF] and GCTA [START_REF] Yang | GCTA: a tool for genome-wide complex trait analysis[END_REF]. Other programs allow the detection of HBD segments in individuals such as PLINK runs of homozygosity (ROHs), BEAGLE [START_REF] Browning | High-resolution detection of identity by descent in unrelated individuals[END_REF] and IBDLD [START_REF] Han | Using identity by descent estimation with dense genotype data to detect positive selection[END_REF]. None of these applications however provide an integrative solution to exploit inbreeding information.

FSuite is a user friendly pipeline composed of several functions integrating FEstim software [START_REF] Leutenegger | Estimation of the inbreeding coefficient through use of genomic data[END_REF]. It estimates f and proposes population genetic statistics that are not available in other software, such as detecting inbred individuals, and inferring parental mating types. A homozygosity mapping statistic with heterogeneity, HFLOD, is proposed, that gives more importance to HBD segments on individuals with small f [START_REF] Leutenegger | Using genomic inbreeding coefficient estimates for homozygosity mapping of rare recessive traits: application to Taybi-Linder syndrome[END_REF]. With these features, it is easy to perform the HBD-GWAS strategy [START_REF] Genin | Could Inbred Cases Identified in GWAS Data Succeed in Detecting Rare Recessive Variants Where Affected Sib-Pairs Have Failed?[END_REF], i.e. homozygosity mapping on inbred cases from genome-wide association study (GWAS), in order to detect Mendelian sub-entities of complex diseases. FSuite also provides graphical outputs to facilitate interpretations of homozygosity mapping results.

METHODS

FSuite needs files in PLINK format: a "map" and a "ped" file, or PLINK binary files. They should contain genome-wide data for the 22 autosomes. If the sample is large enough, FSuite can estimate allele frequencies. Otherwise, frequencies estimated on a reference sample should be furnished in a PLINK "frq" format.

Creation of random submaps

FEstim uses a hidden Markov model (HMM) to model the HBD process of an individual. It requires the markers to be in minimal linkage disequilibrium (LD), which is not the case of actual dense genetic data. Running FEstim on multiple random sparse genome maps (submaps) has been proposed to remove LD [START_REF] Leutenegger | Consanguinity around the world: what do the genomic data of the HGDP-CEPH diversity panel tell us?[END_REF]. FSuite allows the creation of such submaps, based on genetic or physical positions. The default option is to create 100 random submaps from the map file, with one marker every 0.5 cM. An alternative option is to randomly select a marker between recombination hotspots [START_REF] Mcvean | The fine-scale structure of recombination rate variation in the human genome[END_REF][START_REF] Winckler | Comparison of fine-scale recombination rates in humans and chimpanzees[END_REF].

Population genetic studies

The FEstim model depends on 2 parameters f and a, where f is the individual inbreeding coefficient and 1/(a(1-f)) is the expected length of HBD segments (here cM), that are estimated by maximum likelihood. When several submaps are considered, FSuite estimates the inbreeding coefficient as the median value of the f estimates obtained on the different maps. FSuite also fixes FEstim HMM parameters to compute the likelihood of different mating types. These likelihoods can be used for: 1) inferring an individual as inbred by comparing the maximized likelihood with the one to be outbred with a likelihood ratio test, 2) estimating the proportion of mating types in a population, and 3) estimating the probability for an individual to be offspring of different mating types [START_REF] Leutenegger | Consanguinity around the world: what do the genomic data of the HGDP-CEPH diversity panel tell us?[END_REF]. The ones considered in FSuite are first cousin, second cousin, double first cousin, avuncular and unrelated. FSuite outputs the median p-values/probabilities obtained on the multiple submaps.

Homozygosity mapping and HBD-GWAS strategy

HBD posterior probabilities and homozygosity mapping score (FLOD) are calculated per individual at each marker. If a marker is present in several submaps, then an average of the HBD posterior probabilities and FLOD are calculated. Note that FSuite accepts single individuals and nuclear families. Finally, to allow for heterogeneity at a locus, a heterogeneity FLOD (HFLOD) is maximized over a grid of α values (the proportion of cases linked to this locus).

By default, FSuite performs these steps on cases previously inferred as inbred, to remove non-informative outbred individuals for homozygosity mapping. Applied on GWAS data, this allows performing directly HBD-GWAS strategy [START_REF] Genin | Could Inbred Cases Identified in GWAS Data Succeed in Detecting Rare Recessive Variants Where Affected Sib-Pairs Have Failed?[END_REF].

Graphical outputs

For a quicker and easier understanding of the results, FSuite produces some graphical outputs generated with R or Circos [START_REF] Krzywinski | Circos: an information aesthetic for comparative genomics[END_REF], in addition to text output files (see Figure 1 for different examples). Note that graphs similar to the ones of Figures 1-B and1-C can also be plotted for runs of homozygosity detected by PLINK.

Application to whole exome sequencing data

Using submaps has shown good results on SNP data. Whole exome sequencing data are now more and more available, but do not uniformly cover the genome. However, we observed through simulation studies, that FSuite also provides accurate f estimates when applied to SNP genotypes extracted from exome data (see supplementary information).

Requirement and documentation

FSuite needs Perl, R and some of its packages, Circos (optional), PLINK and Merlin [START_REF] Abecasis | Merlin--rapid analysis of dense genetic maps using sparse gene flow trees[END_REF] to be installed on your computer. FSuite pipeline includes a detailed documentation and a simulated example dataset. It also includes an example dataset of 5 cases affected with a rare disease genotyped on 180,160 SNPs, among whom 4 are inbred.

DISCUSSION

FSuite is a user friendly pipeline developed for population genetic studies, rare disease studies and multifactorial disease studies. It implements a set of statistics that are not available in other existing genetic software. It can be used on SNP data but also on whole exome sequencing data. FSuite does not model LD, but takes it into account by using several sparse submaps. It provides robust results, and can be applied on small datasets (where modeling LD is not accurate) if reference allele frequencies are available. In addition, on larger datasets, it has been shown that this strategy gives less bias than available single-points estimates, ROH based estimates, and HMMs modeling LD (Gazal, et al., in press). FSuite offers the possibility to exploit inbreeding information derived from genomic data and will help investigators to better explore their SNP-chip and exome data.

Figure 1 .

 1 Figure 1. FSuite graphical outputs. These graphs were generated with FSuite example dataset. Graph (A) shows the genome-wide HFLOD plot. Dots are HFLOD at each marker. The solid line is a moving average, to remove the impact of a submap with a false positive signal. In graph (B), the HBD segments are represented genome-wide as a Circos plot. The arrow pinpoints chromosome 20, where the maximum HFLOD is reached. Graph (C) shows the HBD segments for each case on this chromosome.

ACKNOWLEDGEMENTS

SG is funded by the plateforme de génomique constitutionnelle (Faculté de médecine, Univ Paris-Diderot, Paris, France).

Resources

1000 Genomes haplotypes and frequencies were downloaded at http://mathgen.stats.ox.ac.uk/impute/data_download_1000G_phase1_integrated.html. EVS frequencies were downloaded at http://evs.gs.washington.edu/EVS/.

Supplementary information: FSuite accuracy on whole exome sequencing data

Accuracy of FSuite on whole exome sequencing (WES) data was checked by simulating 100 offspring of first-cousin using European (EUR) 1000 Genomes (1000G) haplotypes.

Simulations

Individual genomes of 1,000 first-cousin offspring were simulated by gene-dropping with Genedrop program of MORGAN2.9 (www.stat.washington.edu/thompson/Genepi/MORGAN).

To have realistic genome patterns, we downloaded 758 available EUR 1000G autosomal haplotypes phased with SHAPEIT version 2 [START_REF] Delaneau | Improved whole-chromosome phasing for disease and population genetic studies[END_REF] as reference haplotypes.

They were randomly drawn without replacement for each chromosome and were assigned to founders to construct the genotype data of each individual of the family.

We kept 1000G variants coming from two map designs. First, we kept the variants present in the Affymetrix 250K chip. Second, to simulate a WES dataset, we kept variants that are referenced in the exome variant server database (EVS), and that have a minor allele frequency (MAF) ≥ 5% in EUR of 1000G to keep frequent polymorphisms. This cutoff is motivated by the fact that variants with small MAF are numerous but only informative for very few individuals. So, if these variants are selected in a submap, they will bring many frequent uninformative homozygous genotypes. In addition, the frequency of rare alleles is difficult to estimate. A total of 316,963 variants from 1000G were thus kept for this simulation study: 248,290 present in Affymetrix 250K chip, 71,206 in EVS and 2,533 common to both.

To define true HBD in the simulated inbred individuals, founder haplotype labels of the 316,963 variants were used. The true inbreeding coefficient (f true ) of each of the 1,000 first-cousin offspring was calculated by dividing the genome length that is HBD (in cM) by the total genome length obtained by adding the genetic distance between the first and the last variants on each autosome (in cM).

Methods

Different submap selections were performed on data using Affymetrix 250K markers, or our selection of frequent WES polymorphisms. Three methods were compared:

1) FSuite_dft: FSuite with default option (creating 100 random submaps with one marker every 0.5 cM)

2) FSuite_1hot: FSuite with 1 random submap created from recombination hotspots.

3) FSuite_100hot: FSuite with 100 random submaps created from recombination hotspots.

All these methods used the allele frequencies of EUR of 1000G.

Simulation results

For each method, we plotted the f estimates obtained using either Affymetrix 250K data or WES data as a function of the f true values for the 1,000 first-cousin offspring (Figure S1). The results obtained with Affymetrix 250K data and WES data were very similar, whatever the submap selection method, suggesting that FSuite can be used with WES data that do not uniformly cover the genome.

Note that the best correlations were obtained for the multiple submaps based on recombination hotspots, followed by the default option and finally the single submap based on recombination hotspots.


Aymeric Pellé
email: aymeric.pelle@inria.fr

Monique Teillaud
email: monique.teillaud@inria.fr

Periodic meshes for the CGAL library

Keywords: periodic mesh, flat torus, Delaunay refinement, software

This work is motivated by the need for a 3D periodic mesh generator in various fields, such as material engineering or modeling of nano-structures. In this paper, we present a software package answering this need, and which will be made publicly available in the open source CGAL library.

Introduction

CGAL, the Computational Geometry Algorithms Library ( ), is an open source software library offering robust and efficient implementations of various geometric algorithms. It has a number of industrial and academic users in a wide range of fields, like molecular biology, medical modeling, astrophysics, particle physics, geographic information systems, etc.

The library is a collection of C++ software packages with an homogeneous API. Its sticks to the generic programming paradigm through templates, with a style close to the STL programming style [START_REF] Fogel | Generic programming and the CGAL library[END_REF]. Robustness is ensured by following the exact geometric computation paradigm, in which decisions are made on the basis of predicate evaluations that are guaranteed to be exact. Combining exact predicates with filtering techniques allows to gather efficiency and exactness.

CGAL offers in particular a package for volume mesh generation in R 3 , based on Delaunay refinement [START_REF] Jamin | CGALmesh: a Generic Framework for Delaunay Mesh Generation[END_REF][START_REF] Alliez | 3D mesh generation[END_REF]; it uses the 3D triangulation package [START_REF] Pion | 3D triangulations[END_REF]. While that package fulfills the needs of users for many applications, periodic meshes are needed for other fields, such as material engineering and modeling of nano-structures, as confirmed during several multidisciplinary meetings [START_REF]CGAL prospective workshop on geometric computing in periodic spaces[END_REF][START_REF]Subdivide and tile: Triangulating spaces for understanding the world[END_REF][START_REF]Workshop on geometric structures with symmetry and periodicity[END_REF]. To the best of our knowledge, no publicly accessible software is answering that need.

A package for computing periodic Delaunay triangulations in R 3 was introduced in 2009 in CGAL 3.5 [START_REF] Caroli | 3D periodic triangulations[END_REF]. The software was demonstrated by a video [START_REF] Caroli | Video: On the computation of 3D periodic triangulations[END_REF] and has already been used by researchers in other research fields, e.g. condensed matter physics [START_REF] Ogarko | Hierarchical multi-scale modeling[END_REF] and astrophysics [START_REF] Sousbie | The persistent cosmic web and its filament structure I: Theory and implementation[END_REF][START_REF] Van De Weygaert | Alpha, Betti and the megaparsec universe: on the homology and topology of the cosmic web[END_REF].

After quickly reviewing the two CGAL 3D mesh generation and 3D periodic triangulation packages, we present a new software for 3D periodic mesh generation. We show a sample of results it has produced, both on mathematical surfaces and on real world data. The software will be submitted to the CGAL Editorial Board for integration into a future release of the library, under a GPL licence.

CGAL 3D mesh generation

The mesh generation package of CGAL [START_REF] Alliez | 3D mesh generation[END_REF] allows users to generate isotropic tetrahedral meshes of 3D bounded domains. Fig. 1 schematizes its design and the main steps of the computation.

The user defines the domain through an oracle, providing the mesh generator with predicates and constructions that will allow it to query the domain. The domain may consist of multiple components, it may also be subdivided in several subdomains. The behavior of the mesh generator can be customized to meet user-defined requirements, for instance in terms of sizing field or quality (criteria on size or angles of mesh elements). The meshing algorithm is based on Delaunay refinement: in a nutshell, facets in a first step, then cells in a second step, are refined until the criteria are satisfied. The mesh is finally optimized in a later step. All details can be found in [START_REF] Jamin | CGALmesh: a Generic Framework for Delaunay Mesh Generation[END_REF].

CGAL 3D periodic triangulations

Periodic triangulations in R 3 can also be seen as triangulations the 3D flat torus, defined as the quotient space of R 3 under the action of a group generated by three independent translations. The CGAL 3D periodic triangulation package [START_REF] Caroli | 3D periodic triangulations[END_REF] currently allows to handle Delaunay triangulations in the 3D flat torus R 3 /(Z 3 , +), for which the group is generated by three equal length translations along the coordinate axes; its fundamental domain thus is a cube.

Periodic Delaunay triangulations can be constructed and modified incrementally by inserting and removing vertices. The correctness of the algorithm can be ensured in two ways [START_REF] Caroli | Computing 3d periodic triangulations[END_REF]. The default consists in computing in the 27-sheeted covering space R 3 /((3Z) 3 , +), i.e. with 27 copies of each input points, as long as the edges of the triangulation are longer than some value; when all edges are short enough, then the triangulation can safely be converted to R 3 /(Z 3 , +), and all remaining points are inserted without resorting to copies. It was shown that only the first few hundred points are copied, so when the input data set is large, this has a very low impact on the running time. Optionally, the user can choose the heuristic way: it consists in initializing the triangulation with a precomputed triangulation of "dummy" points, chosen in a way to ensure that the above-mentioned condition on edge lengths is always satisfied; then all input points can safely be inserted; the algorithm tries to remove the dummy points at the end of the construction, which is possible if the input point set is sufficiently large and well distributed. Both methods manage to use copies of points only when actually needed, thus allowing for both correctness and efficiency.

Periodic mesh in 3D

Our software uses the above two CGAL packages. As the 3D mesh generation package was designed on top of the CGAL 3D triangulation package, several obstacles must be overcome to use a CGAL 3D periodic triangulation instead.

Some of the difficulties occur at the level of the interface between the two existing packages, and require to modify the CGAL source code, for instance:

• Since a vertex may be associated with several points, copies on an input point, accessing the precise point in the domain must be done through instead of ! (here denotes the triangulation and denotes a vertex handle).

• For the same reason, the periodic criteria must be given more information to access points. This is achieved by passing the triangulation to the criteria as an additional template parameter. Due to lack of space, we omit a few other technical changes to be done. At a higher level, the semantics of the oracle and of the criteria is modified: they must take periodicity into account. For instance, when evaluating the predicate for a segment that intersects a facet of the fundamental domain, it may be the case that does not intersect the surface in the domain, but its translated copy that intersects the opposite domain facet does intersect the surface in the domain (see Fig. 2 for an illustration on a simple shape in 2D). So, for such segments, if the first call of the predicate does not find an intersection, it must be called a second time with the appropriate translated image of the segment. The heuristic option presented in Sec. 3 is chosen for the underlying periodic triangulation, so no point is duplicated during the construction of the mesh, and only the part of the shape that is included in one fundamental domain is meshed. The periodic mesh generator provides users with the same flexibility as the non-periodic one: it allows to mesh multiple domains and to define shape criteria, as illustrated in Fig. 3 for a triply periodic minimal surface, the Schwarz P surface, containing spheres. (Handling domains with sharp features and optimizations is in progress.)

Fig. 1 .

 1 Fig. 1. The CGAL 3D mesh generator (picture taken from [2])

  s

Fig. 2 .Fig. 3 .

 23 Fig. 2. Periodic oracle

Fig. 4

 4 Fig. 4 to 7 show results obtained on the gyroid, another triply periodic minimal surface. Fig. 8 and Fig. 9 show results obtained on real world examples.

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 456789 Fig. 4. The mesh that is actually computed (only the interior tetrahedra are shown) Fig. 5. Visualization of eight copies of the surface facets of the mesh

Acknowledgments

The authors wish to thank Mikhail Bogdanov for his work on a first prototype implementation.


Pietro Abate
email: pietro.abate@pps.univ-paris-diderot.fr

Roberto Di Cosmo
email: roberto@dicosmo.org

Louis Gesbert
email: ocamlprolouis.gesbert@ocamlpro.com

Fabrice Le Fessant
email: fabrice.lefessant@inria.fr

Stefano Zacchiroli

Using Preferences to Tame your Package Manager

Determining whether some components can be installed on a system is a complex problem: not only it is NP-complete in the worst case, but there can also be exponentially many solutions to it. Ordinary package managers use ad-hoc heuristics to solve this installation problem and choose a particular solution, making extremely difficult to change or sidestep these heuristics when the result is not the one we expect.

When software repositories become complex enough, one gets vastly superior results by delegating dependency handling to a specialised solver, and use optimisation functions (or preferences) to control the class of solutions that are found.

The opam package manager relies on the CUDF pivot format, which allows OCaml users that have a CUDF-compliant solver on their machine to reap the benefits of preferencesbased dependency resolution. Thanks to the solver farm provided by Irill, these benefits are now extended to the OCaml community at large.

In this talk we will present the preferences language and explain how to use it.

Introduction

When installing a new software component or upgrading a previously installed one, users often make implicit assumptions regarding the outcome of such operations. To name one, we generally want to install the latest and greatest available version of any software component. Unfortunately, it is not always possible to do so. For example, installing a recent version of a package p may force to downgrade or remove another package q, which fall short of the general desire of only having the latest and greatest components.

What should a package manager do in the above situation? Remove q to get the latest p, or keep q and get the most recent available version of p which is compatible with q The answer really depends on fine-grained user expectations, and different users may have very different points of view and needs. Somebody might want to install the latest version of p because of a new feature that she needs, while someone else might want to just install any version of p and to keep a recent version of q.

Most package managers use ad-hoc heuristics to handle these cases, sometimes leading to surprising results. Moreover, since these heuristics are often hard coded in the dependency solving algorithms, it is usually quite difficult to adapt them to specific needs.

Following the guidelines of [START_REF] Abate | A modular package manager architecture[END_REF], this problem is handled in an original way in the design of the opam package manager, which is specifically tailored to install, upgrade and remove OCaml libraries and tools. Leveraging 10 years of research in the realm of dependency solving, opam uses the de facto standard format CUDF1 to encode upgrade problems, which is recognised by several specialised dependency solvers, and allows users to express complex scenarios in a well defined and compact way by means of preferences that are passed over to the external solver. To make external solvers easily available to all users, independently of their platform or processing power, we have made available at Irill a solver farm that can be used in the Cloud, thus extending the benefits of preferences based dependency solver to all the OCaml community.

Using preferences to guide component installation

A flexible and expressive means of tailoring the behaviour of a dependency solver is offered by the preference language developed by the Mancoosi project, which combines package selectors and measurement functions to give back to the user a fine control on how dependency solving is performed by the package manager. Here we describe the latest version of the preference language, stabilised during the MISC competition in 2012 [START_REF] Abate | Dependency solving: a separate concern in component evolution management[END_REF].

Package selectors

Package selectors are used to identify specific sets of packages, and denoted by the constant expressions: solution, new, removed, changed, up and down .

Given an initial installation I and a proposed solution S, each of these constant expressions denotes a set of packages that can be made formally precise. For example, solution is the set of all packages installed in S, removed is the set of packages that were in I are no longer in S, and so on.

Measurements

A measurement is an integer valued function that can be applied to any of the above package selectors, plus possibly additional arguments. The measurements available in the preference language are the following: count(X), sum(X, f ), unsat recommends(X), notuptodate(X) where X is a set of packages, and f a property of type int or one of its subtypes. These properties may be core properties defined in the CUDF specification, or properties that are declared in the preamble of the CUDF document.

Measurements are evaluated as follows. Note that some of the definitions make use of the set S (the solution) or U (the universe, that is the set of all known packages): count(X) is #(X) that is, the number of packages in X. sum(X,f ) is Σx∈X x.f that is, the sum over all packages in the set X of the value of their field f . notuptodate(X) #{p ∈ X|∃q ∈ U : p.name = q.name ∧ p.version < q.version} that is, notuptodate(X) is the number of packages in X that are not in the latest known version.

unsat recommends(X) counts the number of disjunctions appearing in the Recommends (a.k.a. "soft/optional dependencies") field of packages in X that are not satisfied by S.

Different versions of this preference language are supported by several CUDF-compatible solvers, and can be directly exploited by opam users to tailor the behaviour of opam to their personalised needs.

see http://www.mancoosi.org/cudf/ and the CUDF format specification[START_REF] Treinen | Description of the cudf format[END_REF]

http://www.cs.uni-potsdam.de/wv/aspcud/

Preferences A user preference is a sequence of measurements of package selectors, each prefixed by either the plus or minus sign, to indicate to the solver whether we desire to see the corresponding measure maximised or minimised in the solution.

For example, an opam user interested in installing a new package at its latest version, but touching as little as possible all other installed packages, may write the following preference expression:

-notuptodate(new),-count(changed)

Such expression can be passed to opam directly on the command line, using --criteria.

Dependency solving in the cloud

To be able to leverage the best dependency solvers, it is important to keep the solver as a modular component, separate from the rest of the package manager, as explained in [START_REF] Abate | Dependency solving: a separate concern in component evolution management[END_REF][START_REF] Abate | A modular package manager architecture[END_REF]. This may lead to some added complexity in setting up a multi-platform package manager like opam, as it is necessary to make sure that the existing dependency solvers (among which aspcud 2 is the most advanced) are available, with an aligned version, on all platforms.

For this reasons, a task force at Irill has been organised, with help from OCamlPro, and in particular by Pierre Chambart, Benjamin Canou and Gregoire Henry, to set up a solver farm which is online at http://cudf-solvers.irill.org.

The solver farm can be easily accessed over the network by opam to perform dependency solving in the Cloud. This can be particularly useful for developers working on devices with limited resources, like the Raspberry Pi. The testing phase has shown that the communication and solving time amount to less than two seconds, even on legacy internet connections. Hence the benefits of personalised dependency solving are now available to all the OCaml community.

Conclusion

In this talk we will present the basic underpinnings of external dependency solving for a package manager, with a special focus on the preferences language and its expressive power; we will also show how to use the CUDF solver farm in conjunction with the opam package manager.

Software availability

The dependency solving toolchain covered by this talk is made of several building blocks, all of which are entirely written in OCaml and available as Free Software:

• opam -the new OCaml package manager, can be found at http://opam.ocamlpro.com/

• DOSE -the comprehensive library for dependency solving and analysis, collecting the works developed during the EDOS and MANCOOSI project, can be found at https: //gforge.inria.fr/projects/dose/ • libCUDF -the reference implementation for manipulating and checking CUDF files, is available at http://www.mancoosi.org/cudf/


ISAE Pierre Roux

Innocuous Double Rounding of Basic Arithmetic Operations

Double rounding occurs when a floating-point value is first rounded to an intermediate precision before being rounded to a final precision. The result of two such consecutive roundings can differ from the result obtained when directly rounding to the final precision. Double rounding practically happens, for instance, when implementing the IEEE754 binary32 format with an arithmetic unit performing operations only in the larger binary64 format, such as done in the PowerPC or x87 floating-point units. It belongs to the folklore in the floating-point arithmetic community that double rounding is innocuous for the basic arithmetic operations (addition, division, multiplication, and square root) as soon as the final precision is about twice larger than the intermediate one. This paper adresses the formal proof of this fact considering underflow cases and its extension to radices other than two.

INTRODUCTION

Floating-point numbers are commonly used to efficiently perform numerical computations, which are then performed with a bounded precision. That is, only a finite number of bits are kept after each arithmetic operation. Multiple choices are available for this precision, typical examples being 24 bits for the IEEE754 binary32 format and 53 for the binary64 format.

Double rounding occurs when a value is first rounded to an intermediate precision before being rounded to a final precision. The result of two such consecutive roundings can differ from the result obtained when directly rounding to the final precision. In radix 10, for instance, 1.495 could be rounded to 1.50 then to 2 if first rounding to a 3-digit precision then to one digit whereas a direct rounding to the nearest value with one digit would give 1. Double rounding practically happens, for instance, when implementing the IEEE754 binary32 [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF] format with an arithmetic unit only performing operations in the larger binary64 format, such as done in the PowerPC or x87 floating-point units.

It belongs to the folklore in the floating-point arithmetic community that double rounding is innocuous for the basic arithmetic operations (addition, subtraction, multiplication, division and square root) as soon as the final precision is about twice larger than the intermediate one. The IEEE754 binary32 and binary64 formats [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF] fulfill this condition. Figueroa [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] published pen-and-paper proofs of these results but only considering radix 2 and completely ignoring underflows. We try to give more general results in this paper. Moreover, such proofs tend to involve quite a number of rather subtle corner cases which make them particularly error prone. We therefore use Coq [START_REF] Bertot | Interactive theorem proving and program development : Coq'Art : the calculus of inductive constructions[END_REF][START_REF]The Coq proof assistant reference manual[END_REF] to ensure no such corner case has been missed. The impact of double rounding on some particular floating-point al-

• Pierre Roux gorithms such as the Veltkamp/Dekker's algorithms has already been studied using Coq [START_REF] Martin-Dorel | Some issues related to double rounding[END_REF], but its innocuousness, under some conditions, on basic arithmetic operations is only briefly mentioned.

More formally, given a floating-point format F, a rounding • 1 in this format and a second rounding • 2 in a larger-precision format, we want to ensure:

∀x, y ∈ F, • 1 (• 2 (x ⋄ y)) = • 1 (x ⋄ y)
under some conditions, for ⋄ ∈ {+, -, ×, /, √ }. That is, the second rounding • 2 is innocuous for basic arithmetic operations. This paper aims at studying under which conditions the previous holds and to formally prove such results using a proof assistant (Coq). It is worth noting that the results are highly dependent on the operation ⋄ and that they may not hold for operations other than {+, -, ×, /, √ }.

Since all the proofs presented in the remaining of the paper have been formally verified within Coq, we only present proof sketches. Although not all subcases are precisely demonstrated, enough is shown to give a good idea of what is going on. The case of directed roundings (toward -∞ or +∞ or roundings to zero or away) is rather simple: the above result holds if and only if • 1 and • 2 are rounding in the same direction. The remaining will therefore focus on the case where • 1 and • 2 are roundings to nearest (possibly with different tie-break rules). We can nevertheless notice that the results still hold (under possibly weaker hypotheses) if • 1 is a directed rounding and • 2 a rounding to nearest. The converse case (• 1 to nearest and • 2 directed) however usually requires slightly stronger hypothesis.

Finally, we do not deal with overflows since they are harmless as soon as any number which can be represented in the smaller precision format F can also be represented in the larger precision, which is a reasonable assumption on exponent ranges.

The remainder of this section introduces the floating-point formalisms used in our proofs. Then, Sections 2, 3, 4, and 5 illustrate our proofs about respectively the multiplication, the addition/subtraction, the square root and the division. For the sake of clarity, the proofs are presented on a formalism not modeling underflows, although our actual Coq proofs handle them and precise requirements on underflow will be given. Eventually, Section 6 summarizes the results and concludes.

All our Coq developments are part of the Flocq library [BM11] available1 at http://flocq.gforge.inria.fr/.

Floating-point Formats

We first formally define some floating-point formats used throughout the remaining of the paper. A generic format is first defined, then particular cases: a fixed precision format without underflow, a format with gradual underflows and one with abrupt underflows. All these formats are formally defined in the Coq library for floating-point arithmetic Flocq [START_REF] Boldo | Flocq: A Unified Library for Proving Floating-point Algorithms in Coq[END_REF].

For this purpose, we need a discrete logarithm.

• 3

Definition 1. Given β ∈ Z, β ≥ 2, we define the discrete logarithm ln β (x) for radix β of a real number x ∈ R \ {0} as the unique value e x in Z satisfying

β ex-1 ≤ |x| < β ex .
A generic format F ϕ is defined by a function ϕ : Z → Z giving, for all numbers sharing the same discrete logarithm, the precision they must be encoded with.

Definition 2. Given ϕ : Z → Z, a value x ∈ R is said to be in the generic format

F ϕ if x = 0 or if there is some m ∈ Z such that x = m β ϕ(ln β (x)) .
m is then called the mantissa of x and ϕ(ln β (x)) its canonical exponent. In the Flocq library, the proposition generic format β ϕ x means that x ∈ F ϕ . This generic format then allows us to define more concrete formats. Given some p ∈ Z, p ≥ 1, the format FLX p models precision p floating-point numbers without underflow.

Definition 3. Given p ∈ Z, p ≥ 1, the format FLX p is defined as F ϕ with ϕ : e → e -p (denoted FLX exp p in Flocq).

Although this format remains mostly theoretical since hardware floating-point numbers can underflow, it is still a good model of what happens with actual floatingpoint values when no underflow occurs.

Given p ∈ Z, p ≥ 1, and e min ∈ Z, the format FLT p,emin models precision p floating-point numbers with underflow handled by subnormals of exponent e min . Definition 4. Given p, e min ∈ Z, p ≥ 1, the format FLT p,emin is defined as F ϕ with ϕ : e → max(e -p, e min ) (denoted FLT exp e min p in Flocq).

Remark 5. FLT 24,-149 with radix β = 2 corresponds to the IEEE754 binary32 format and FLT 53,-1074 to binary64 (neglecting overflows).

Finally, given p ∈ Z and e min ∈ Z, the format FTZ p,emin models precision p floating-point numbers with underflow handled by flushing either to zero or to the smallest normal number. Definition 6. Given p, e min ∈ Z, the format FTZ p,emin is defined as F ϕ with ϕ : e → e -p when e ≥ e min + p and e → e min + p -1 otherwise (denoted FTZ exp e min p in Flocq).

The ϕ functions defining FLX, FLT and FTZ formats are illustrated on Figure 1. It only remains to define roundings.

Definition 7. A rounding on a floating point format F is a function • : R → F monotone and equal to the identity on F, that is

∀x, y ∈ R, x ≤ y ⇒ •(x) ≤ •(y) ∀f ∈ F, •(f ) = f.
This corresponds to the hypothesis Valid rnd rnd in Flocq while the rounding function • on F ϕ is denoted by round β ϕ rnd. Definition 8. A rounding to nearest is a rounding • that returns a nearest value in F:

∀x ∈ R, ∀f ∈ F, |x -•(x)| ≤ |x -f | .
They are of the form round β ϕ (Znearest ) in Flocq.

Remark 9. In case x lies halfway between two consecutive values in F, •(x) is not uniquely defined as it can be any of these two consecutive values. This choice is called a tie-break rule. Since most of our proofs are valid regardless of the tie-break rule (i.e., for any rounding to nearest) we don't extend on this point.

MULTIPLICATION

Proofs for the addition/subtraction being rather involved, let us begin with the multiplication for which proofs are much simpler.

General Case

Figueroa [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] proved double rounding to be innocuous for the multiplication if it is performed with a precision at least twice as large. Although his proof only addressed radix 2, this holds for any radix.

Theorem 10. ( [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF], double round mult FLX in our Coq development) For p 1 , p 2 ∈ Z, if p 2 ≥ 2p 1 , then for • 1 and • 2 any roundings (for instance, toward -∞ or to nearest with any tie) respectively in FLX p1 and FLX p2 :

∀x, y ∈ FLX p1 , • 1 (• 2 (x × y)) = • 1 (x × y) .
This in particular applies to binary32/64 when no underflow occurs.

Remark 11. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] This bound 2p 1 is optimal as shown by the following counterexample. With radix β = 2, for p 1 = 4, p 2 = 7, we have for x = y = 1.101 2

• 1 (• 2 (x × y)) = • 1 (• 2 (10.101001 2 )) = • 1 (10.10100 2 ) = 10
.10 2 which differs from • 1 (10.101001 2 ) = 10.11 2 for • 1 and • 2 appropriate roundings to nearest (for instance, the common rounding to nearest with tie-break to even).

• 5

Proof. x and y having p 1 significand digits, x × y has at most 2p 1 ≤ p 2 digits hence

• 2 (x × y) = x × y.
The previous result still applies to binary32/64 even in case of underflow. Indeed, proofs for the FLX, FLT and FTZ formats (double round mult {FLX,FLT,FTZ} in our Coq development) are just corollaries of a single proof (double round mult) carried out on the generic format F ϕ . This general proof shows innocuity of double rounding of multiplication under the assumption ∀e x , e y ∈ Z, ϕ 2 (e x + e y ) ≤ ϕ 1 (e x ) + ϕ 1 (e y ) ∧ ϕ 2 (e x + e y -1) ≤ ϕ 1 (e x ) + ϕ 1 (e y ), the main argument being the same than in the above proof.

Odd Radix

It is worth noting that a much better result can be obtained for odd radices β, since the hypothesis p 2 ≥ 2p 1 in Theorem 10 can be replaced with the weaker p 2 ≥ p 1 when β is odd. This nice result is, at least currently2 , certainly perfectly useless. We nevertheless chose to develop it for the sake of exhaustivity and considering it did not implied a huge overhead on our Coq development.

Theorem 12. (double round mult beta odd FLX) When β is odd, for p 1 , p 2 ∈ Z, if p 2 ≥ p 1 , then for • 1 and • 2 roundings to nearest, with any tie, respectively in

FLX p1 and FLX p2 ∀x, y ∈ FLX p1 , • 1 (• 2 (x × y)) = • 1 (x × y) .
To give an idea of the proof, we first need a few additional definitions and lemmas.

Definition 13. • ↓ ϕ denotes rounding towards -∞, i.e., • ↓ ϕ (x) is the largest floatingpoint value in F ϕ below x:

• ↓ ϕ : x ∈ R → max {f ∈ F ϕ | f ≤ x}. Definition 14. ulp ϕ (x) denotes the unit in last place of x in format F ϕ : ulp ϕ (x) := β ϕ(ln β (x)) .
Definition 15. midp ϕ (x) denotes the following value, called midpoint

midp ϕ (x) := • ↓ ϕ (x) + ulp ϕ (x)
2 as it lies halfway between two consecutive floating-point values around x.

Lemma 16. (double round lt mid further place) Given two floating-point formats F ϕ1 and F ϕ2 , two roundings to nearest, with any tie, • 1 and • 2 in these formats, for all

x ∈ R, if x > 0, ϕ 2 (ln β (x)) < ϕ 1 (ln β (x)), ϕ 1 (ln β (x)) ≤ ln β (x) and x < midp ϕ1 (x) - ulp ϕ2 (x) 2 then • 1 (• 2 (x)) = • 1 (x) . • Pierre Roux x • ↓ ϕ 1 (x) midp ϕ 1 (x) ulp ϕ 1 (x) 2 ulp ϕ 2 (x) 2 • 2 (x) Fig. 2: Illustration of Lemma 16. Since x is below midp ϕ 1 (x), • 1 (x) = • ↓ ϕ 1 (x) and since x is below midp ϕ 1 (x) -ulp ϕ 2 (x) /2 then • 2 (x) is also below midp ϕ 1 (x), implying • 1 (• 2 (x)) = • ↓ ϕ 1 (x). x 0 x 1 x 2 x 3 m ulp ϕ 1 (x 0 ) 2 ulp ϕ 2 (x 1 ) 2
Fig. 3: Illustration of Lemma 17. In radix β = 3, with ϕ 2 = ϕ 1 -1, x 0 and x 3 are consecutive values in Fϕ 1 while x 0 , x 1 , x 2 and x 3 are consecutive values in Fϕ 2 . The midpoint m in Fϕ 1 is also a midpoint in

Fϕ 2 : m = x 0 + ulp ϕ 1 (x 0 ) 2 = x 1 + ulp ϕ 2 (x 1 ) 2
. This property relies on β being odd.

This is illustrated on Figure 2.

Proof. Since ϕ 2 (ln β (x)) ≤ ϕ 1 (ln β (x)), we get • 2 (x)-• ↓ ϕ1 (x) = • 2 (x) -• ↓ ϕ1 (x) . Moreover, • 2 being a rounding to nearest, we get for all y ∈ R, |• 2 (y) -y| ≤ ulp ϕ 2 (y) 2
and the hypothesis x < midp ϕ1 (x)-

ulp ϕ 2 (x) 2 gives us x -• ↓ ϕ1 (x) < ulp ϕ 1 (x)-ulp ϕ 2 (x) 2 . Thus • 2 (x) -• ↓ ϕ1 (x) ≤ |• 2 (x) -x| + x -• ↓ ϕ1 (x) < ulp ϕ 1 (x) 2
.

• 1 being a rounding to nearest, • 1 (• 2 (x)) and • 1 (x) are then both equal to • ↓ ϕ1 (x). Contrary to the previous lemmas which were valid for any radix β, the following result shows a particularly interesting property of odd radices.

Lemma 17. (midpoint beta odd remains) When β is odd, for all ϕ 1 , ϕ 2 : Z → Z, for all x ∈ F ϕ1 , if ϕ 2 (ln β (x)) ≤ ϕ 1 (ln β (x)), then there exists an x ′ ∈ F ϕ2 such that x +

ulp ϕ 1 (x) 2 = x ′ + ulp ϕ 2 (x) 2
.

This result is illustrated on Figure 3. Intuitively, this means that, with odd radices, if a real number x +

ulp ϕ 1 (x) 2
is a midpoint for some precision, it is also a midpoint in any larger precision. This is proved by induction on ϕ 1 (ln β (x))ϕ 2 (ln β (x)) ∈ N and makes it possible to prove the following lemma.

Lemma 18. (neq midpoint beta odd) When β is odd, for all ϕ 1 , ϕ 2 : Z → Z, for all x ∈ R, if ϕ 2 (ln β (x)) ≤ ϕ 1 (ln β (x)) and x = midp ϕ1 (x) then • 1 (• 2 (x)) = • 1 (x)
for any rounding to nearest, with any tie, • 1 and • 2 respectively in F ϕ1 and F ϕ2 .

This means that with odd radices, double rounding is innocuous for all real values which are not midpoints.

• 7

Proof. Lemma 16 concludes when x < midp ϕ1 (x)-

ulp ϕ 2 (x) 2
. Let us then assume

midp ϕ1 (x) - ulp ϕ2 (x) 2 ≤ x < midp ϕ1 (x) .
From Lemma 17, there exist an x ′ ∈ F ϕ2 such that x ′ = midp ϕ1 (x) -

ulp ϕ 2 (x) 2 and x ′ ≤ x < x ′ + ulp ϕ2 (x) 2 hence • 2 (x) = x ′ . Since • ↓ ϕ1 (x) ≤ x ′ < midp ϕ1 (x), then • 1 (x ′ ) = • ↓ ϕ1 (x) = • 1 (x). The proof is similar for x > midp ϕ1 (x).
A last lemma is then needed to prove the theorem.

Lemma 19. (float neq midpoint beta odd) When β is odd, for all x ∈ R, ϕ : Z → Z, if there exist m, e ∈ Z such that x = mβ e , then x = midp ϕ (x)

Proof. If e ≥ ϕ(ln β (x)), then both x and • ↓ ϕ (x) are multiples of β ϕ(ln β (x)) whereas ulp ϕ (x) 2

is not. Otherwise, Lemma 17 allows us to conclude by a similar reasoning.

Proof (Theorem 12). Since x × y = (m x × m y ) β ex+ey , with m x , m y and e x , e y the mantissas and exponents of x and y, Lemmas 18 and 19 make it possible to conclude.

ADDITION

General Case

Figueroa [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] proved double rounding to be innocuous for the addition if it is performed with a precision at least strictly twice as large. Although his proof only addressed radix 2, this holds for any radix.

Theorem 20. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] (double round plus FLX) For p 1 , p 2 ∈ Z, if p 2 ≥ 2p 1 +1, then for • 1 and • 2 roundings to nearest, with any tie, respectively in FLX p1 and FLX p2 :

∀x, y ∈ FLX p1 , • 1 (• 2 (x + y)) = • 1 (x + y) .
This applies to IEEE754 binary32/64 when no underflow occurs. Remark 21. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] This bound 2p 1 + 1 is optimal as shown by the following counterexample. With radix β = 2, for p 1 = 4, p 2 = 8, we have for x = 1.001 2 and y = 0.00001111 2

• 1 (• 2 (x + y)) = • 1 (• 2 (1.00101111 2 )) = • 1 (1.0011000 2 ) = 1.010 2
which differs from • 1 (1.00101111 2 ) = 1.001 2 when • 1 and • 2 are roundings to nearest with tie-break to even.

We need the following lemma to prove the previous theorem.

Lemma 22. For all ϕ : Z → Z, x ∈ F ϕ , y ∈ R, if x > 0, y > 0 and ln β (y) ≤ ϕ(ln β (x)) -2, then 0 < (x + y) -• ↓ ϕ (x + y) < β ϕ(ln β (x))-2 .
• Pierre Roux

Proof. Since ln β (y) ≤ ϕ(ln β (x)) -2, we get y < ulp ϕ (x) = ulp ϕ (x + y), hence

• ↓ ϕ (x + y) = x, that is (x + y) -• ↓
ϕ (x + y) = y which implies the result. Proof (Theorem 20). We will only consider the case x > 0, y > 0 and y ≤ x, other cases being pretty similar or easy to deduce. Either ln β (y) ≥ ϕ 1 (ln β (x))-1 in which case x+y has at most 2p 1 +1 significand digits, which means that • 2 (x + y) = x + y, or ln β (y) ≤ ϕ 1 (ln β (x)) -2 in which case the result follows from Lemmas 16 and 22.

The previous result still applies to IEEE754 binary32/64 even in case of underflow (double round plus FLT), with a similar proof thanks to the generic format F ϕ . Indeed, proofs for formats with or without underflow are derived as immediate corollaries of a proof on the generic format.

Under the assumption β ≥ 3, the bound 2p 1 + 1 can be replaced by 2p 1 in Theorem 20. The proof is similar (double round plus beta ge 3 FLX).

Remark 23. This bound 2p 1 is optimal as shown by the following counterexample. With radix β = 10, for p 1 = 4, p 2 = 7, we have for x = 1.001 10 and y = 0.0004995 10

• 1 (• 2 (x + y)) = • 1 (• 2 (1.0014995 10 )) = • 1 (1.001500 10 ) = 1.002 10
which differs from • 1 (1.0014995 10 ) = 1.001 10 when • 1 and • 2 are roundings to nearest with tie-break to even.

Finally, we can notice that the same results can be immediately deduced for the subtraction.

Odd Radix

As for the multiplication, a way better bound is obtained for odd radices.

Theorem 24. (double round plus beta odd FLX) When β is odd, for p 1 , p 2 ∈ Z, if p 2 ≥ p 1 , then for • 1 and • 2 roundings to nearest, with any tie, respectively in FLX p1 and FLX p2 :

∀x, y ∈ FLX p1 , • 1 (• 2 (x + y)) = • 1 (x + y) .
Proof. Denoting m x , m y and e x , e y the mantissas and exponents of x and y, x+y = mβ e with e = min(e x , e y ) and m = m x β ex-min(ex,ey) +m y β ey-min(ex,ey) ∈ Z. Lemmas 18 and 19 make it possible to conclude.

SQUARE ROOT

General Case

Figueroa [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] proved double rounding to be innocuous for the square root if it is performed with a precision larger than twice the original precision plus two. Although his proof only addressed radix 2, this holds for any radix.

Theorem 25. [Fig95] (double round sqrt FLX) For p 1 , p 2 ∈ Z, if p 2 ≥ 2p 1 +2, then for • 1 and • 2 roundings to nearest respectively in FLX p1 and FLX p2 ∀x ∈ FLX p1 , • 1 • 2 √ x = • 1 √ x .
This applies to IEEE754 binary32/64 when no underflow occurs.

• 9

Remark 26. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] This bound 2p 1 + 2 is optimal as shown by the following counterexample. With radix β = 2, for p 1 = 4, p 2 = 9, we have for x = 0.1111 2

• 1 • 2 √ x = • 1 (• 2 (0.1111011111 . . . 2 )) = • 1 (0.111110000 2 ) = 1.000 2
which differs from • 1 (0.1111011111 . . . 2 ) = 0.1111 2 when • 1 and • 2 are roundings to nearest with tie-break to even.

Proof (Theorem 25). If √ x < midp ϕ1 ( √ x)- ulp ϕ 2 ( √ x)

2

, the result follows from Lemma 16. Similarly, the result holds when

√ x > midp ϕ1 ( √ x) + ulp ϕ 2 ( √ x)

2

. We will now show that x cannot lie between those two bounds. Let us first denote

u 1 := ulp ϕ1 ( √ x), u 2 := ulp ϕ2 ( √ x), a := • ↓ ϕ1 ( √ x), b := u1-u2 2 and b ′ := u1+u2 2 . Then assuming a + b ≤ √ x ≤ a + b ′
we get

a 2 + u 1 a -u 2 a + b 2 ≤ x ≤ a 2 + u 1 a + u 2 a + b ′2 .
a 2 +u 1 a as well as x can be divided by u 2 1 . However -u 2 a+b 2 > 0 and u 2 a+b ′2 < u 2 1 which is absurd.

The previous result still applies to IEEE754 binary32/64 even in case of underflow (double round sqrt FLT).

Under the assumption β ≥ 4, the bound 2p 1 + 2 can be replaced by 2p 1 + 1 in Theorem 25. The proof is similar (double round sqrt beta ge 4 FLX).

Remark 27. This bound 2p 1 + 1 is optimal as shown by the following counterexample. With radix β = 10, for p 1 = 4, p 2 = 8, we have for x = 0.9999 10

• 1 • 2 √ x = • 1 (• 2 (0.999949998 . . . 10 )) =
• 1 (0.99995000 10 ) = 1.000 10 which differs from • 1 (0.999949998 . . . 10 ) = 0.9999 10 when • 1 and • 2 are roundings to nearest with tie-break to even.

Odd Radix

As for the multiplication and addition, a way better bound is obtained for odd radices.

Theorem 28. (double round sqrt beta odd FLX) When β is odd, for p 1 , p 2 ∈ Z, if p 2 ≥ p 1 , then for • 1 and • 2 roundings to nearest respectively in FLX p1 and FLX p2 :

∀x ∈ FLX p1 , • 1 • 2 √ x = • 1 √ x .
Proof. According to Lemma 18, it is enough to prove that √ x = midp ϕ1 ( √ x).

Assuming the contrary, we get

x = r 2 + ru + u 2 4
where r := • ↓ ϕ1 ( √ x) and u := ulp ϕ1 ( √ x). Then, when ϕ 1 (ln β (x)) ≥ 2ϕ 1 (ln β ( √ x)),

we have x, r 2 and ru multiples of u 2 whereas u 2 4 is not. Otherwise, Lemma 17 makes it possible to conclude with a similar reasoning.

• Pierre Roux

DIVISION

Even Radix

Figueroa [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] proved double rounding to be innocuous for the division if it is performed with a precision at least twice as large. Although his proof only addressed radix 2, this holds for any even radix.

Theorem 29. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] (double round div FLX) When β is even, for p 1 , p 2 ∈ Z, if p 2 ≥ 2p 1 , then for • 1 and • 2 roundings to nearest, with any tie, respectively in FLX p1 and FLX p2 :

∀x, y ∈ FLX p1 , y = 0 ⇒ • 1 (• 2 (x/y)) = • 1 (x/y) .
This applies to IEEE754 binary32/64 when no underflow occurs.

Remark 30. [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] This bound 2p 1 is optimal as shown by the following counterexample. With radix β = 2, for p 1 = 4, p 2 = 7, we have for x = 1.000 2 and y = 0.1111 2

• 1 (• 2 (x/y)) = • 1 (• 2 (1.000100010 . . . 2 )) = • 1 (1.0001000 2 ) = 1.000 2
which differs from • 1 (1.000100010 . . . 2 ) = 1.001 2 when • 1 and • 2 are roundings to nearest with tie-break to even.

Proof. Denoting midp ϕ1 (x/y) with m and ulp ϕ 2 (x/y) 2

with u, from Lemma 16, the result holds when x/y < m -u and similarly when x/y > m + u. When m -u ≤ x/y < m, an argument similar to the one used in proof of Theorem 25 allows to conclude and similarly when m < x/y ≤ m + u. Finally, when x/y = m, since β is even, x/y ∈ F ϕ2 as soon as ϕ 2 (x/y) < ϕ 1 (x/y) which is implied by p 2 ≥ 2p 1 and p 1 ≥ 1.

The previous result still applies to IEEE754 binary32/64 even in case of underflow (double round div FLT).

Odd Radix

The result does not hold for odd radices.

Remark 31. The result does not hold for odd radices, even for arbitrarily large p 2 , as shown by the following counterexample. With radix β = 3, for p 1 = 4, p 2 ≥ 4, p 2 odd, for x = 2001 3 and y = 2 3

• 1 (• 2 (x/y)) = • 1 (• 2 (1000.111 . . . 3 )) = • 1 (1000.1 . . . 1 3 ) = 1000 3
which differs from • 1 (1000.111 . . . 3 ) = 1001 3 when • 1 and • 2 are roundings to nearest with tie-break to even.

However, the result can be recovered by restricting the choice of tie-break rules to directed rules such as round to nearest with tie-break away.

Theorem 32. (double round div rna FLX) When β is odd, for p 1 , p 2 ∈ Z, if p 2 ≥ p 1 , then for • A 1 and • A 2 roundings to nearest with tie-break away, respectively in FLX p1 and FLX p2 :

∀x, y ∈ FLX p1 , y = 0 ⇒ • A 1 • A 2 (x/y) = • A 1 (x/y) .
• 11 Table I. Hypotheses under which double rounding is innocuous for the FLX format (β is the radix and p 1 and p 2 are the two precisions involved (c.f., Section 1 for formal definitions)). Division for odd radices is only valid with tie-break away. Proof. According to Lemma 18, the result holds when x/y is not a midpoint. It also holds when it is, since the directed tie-break rule ensures that midpoints are always rounded in the same direction.

β = 2 even β ≥ 4 odd β +, -p 2 ≥ 2p 1 + 1 p 2 ≥ 2p 1 p 2 ≥ p 1 × p 2 ≥ 2p 1 p 2 ≥ 2p 1 p 2 ≥ p 1 / p 2 ≥ 2p 1 p 2 ≥ 2p 1 p 2 ≥ p 1 √ p 2 ≥ 2p 1 + 2 p 2 ≥ 2p 1 + 1 p 2 ≥ p 1

CONCLUSION

Although all the results were illustrated throughout the paper on the FLX floatingpoint format (no underflows), they remain valid under mild conditions on the FLT (gradual underflow) and FTZ (abrupt underflow) formats with underflows. Indeed, all our Coq proofs are carried out at the level of F ϕ generic formats and results for the FLX, FLT and FTZ formats with or without underflows are just corollaries immediatly derived from the former. The precise hypotheses under which double rounding has been proved innocuous are summarized in Tables I, II and III.

Once the proofs, originally expressed in terms of bit patterns [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF], have been translated to the more algebraic form seen throughout the paper, Coq proofs appeared rather easy to perform. A key factor in this easyness was the high genericity of the definitions in the Flocq library [START_REF] Boldo | Flocq: A Unified Library for Proving Floating-point Algorithms in Coq[END_REF]. For instance, most of our proofs are done for any rounding to nearest, for any tie, by just using the definition of rounding to nearest available in Flocq. It should also be noticed that the library appeared pretty complete since no definition and only a few (13) lemmas were added to it.

All our Coq developments (6.9 kloc) are part of the Flocq library [START_REF] Boldo | Flocq: A Unified Library for Proving Floating-point Algorithms in Coq[END_REF] which is available 1 under an open source license at http://flocq.gforge.inria.fr/.

Thus, we have not only formally proved previously-known results on innocuousness of double rounding of basic arithmetic operations [START_REF] Samuel | When is Double Rounding Innocuous?[END_REF] but also extended them to radices other than two (such as radix 10 introduced in the last revision of the IEEE754 norm [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF]) and underflow cases (either gradual or abrupt). To the extent of author knowledge, although they belong to the folklore in the floatingpoint arithmetic community, no proof, neither pen and paper nor mechanical, can be found in the literature. In particular, it is now formally proved that double rounding of addition/subtraction, multiplication, division, and square root, is innocuous for IEEE754 binary32 format when using binary64 as intermediate precision. This could, for instance, be used in the verified C compiler Compcert [START_REF] Boldo | A Formally-Verified C Compiler Supporting Floating-Point Arithmetic[END_REF] to perform constant folding while cross compiling for a binary32 target on a binary64 host. This may also be of interest for some processors, such as the one using the POWER (ancestor of PowerPC) instruction set, which have binary64 arithmetic instructions but lack binary32 ones. The RAD6000, a radiation-hardened board used onboard various spacecrafts, belongs to this category.

Fig. 1 :

 1 Fig.1: Illustration of the functions ϕ defining the formats FLX, FLT and FTZ. The difference (eϕ(e)) between the diagonal (doted line) and the graph of ϕ is the size of the mantissa for numbers of magnitude e in format Fϕ (if non positive, numbers can only be encoded as 0). These functions are the same for normal numbers (e ≥ e min + p) and only differ for subnormal numbers. FLX doesn't model them while FLT has gradual underflow and FTZ abrupt underflow.

Table II .

 II Hypotheses under which double rounding is innocuous for the FLT format, in addition to hypotheses of Table I (e min1 and e min2 are the two minimal exponents involved). Division for odd radices is only valid with tie-break away.

	+, -× / √	β = 2 e min2 ≤ e min1 e min2 ≤ 2e min1 e min2 ≤ e min1 -p 1 -2 (e min2 ≤ e min1 -p 1 -2 ∨ 2e min2 ≤ e min1 -4p 1 -2)	even β ≥ 4 e min2 ≤ e min1 e min2 ≤ 2e min1 e min2 ≤ e min1 -p 1 -2 (e min2 < e min1 -p 1 ∨ 2e min2 ≤ e min1 -4p 1 )	odd β e min2 ≤ e min1 e min2 ≤ e min1 e min2 ≤ e min1 e min2 ≤ e min1

Table III .

 III Hypotheses under which double rounding is innocuous for the FTZ format, in addition to hypotheses of TableI. Division for odd radices is only valid with tie-break away.

		even β	odd β
	+, -	e

min2 + p 2 ≤ e min1 + 1 e min2 + p 2 ≤ e min1 + p 1 × e min2 + p 2 ≤ 2e min1 + p 1 e min2 + p 2 ≤ e min1 + p 1 / e min2 + p 2 < e min1 e min2 + p 2 ≤ e min1 + p 1 √ 2 (e min2 + p 2 ) ≤ e min1 + p 1 ≤ 1 e min2 + p 2 ≤ e min1 + p 1

In files src/appli/Fappli double round.v and examples/Double round beta odd.v, starting with version

2.4.0 of Flocq.Journal of Formalized Reasoning Vol. ?, No. ?, Month Year.

Journal of Formalized Reasoning Vol. ?, No. ?, Month Year.

There have been a couple of computers using a radix

arithmetic, e.g., Setun built in Moscow in the 60s and Ternac developped in New York in the 70s.Journal of Formalized Reasoning Vol. ?, No. ?, Month Year.

ACKNOWLEDGMENTS

The author wants to deeply thank Sylvie Boldo and Guillaume Melquiond for their help for this work and the anonymous reviewers for their comments.