Milan Zdravković

Hervé Panetto

Miroslav Trajanović

SEMANTIC INTEROPERABILITY FOR DYNAMIC PRODUCT-SERVICE BUSINESS ECOSYSTEMS

Dynamic Product-Service Business Ecosystems (DPSBEs) are collaborative forms that aim to handle product-service business actions in a business ecosystem by creating, sharing, using and reusing information across functional teams and enterprise boundaries. DPSBEs depend on the ability of their members' information systems to interoperate. This ability is related to making their semantics explicit and formal, so information can be interpreted and managed both by humans and systems. In this position paper, we establish the scientific background for the semantic interoperability for dynamic product-service business ecosystems.

INTRODUCTION

Dynamic Product-Service Business Ecosystems (DPSBEs) are emerging collaborative forms that aim at organising and handling product-service business actions.

Product-Service System (PSS) is thought [START_REF] Baines | State-of-the-art in product service-systems[END_REF] as "a market proposition that extends the traditional functionality of a product by incorporating additional services. Here the emphasis is on the 'sale of use' rather than the 'sale of product'. The customer pays for using an asset, rather than its purchase, and so benefits from a restructuring of the risks, responsibilities, and costs traditionally associated with ownership".

The PSS concept embraces an integrated view on the products and associated services, supporting networks and infrastructure and design and innovation strategy to achieve and implement this view. It assumes improved competitiveness, better customer satisfaction (more product customization and higher quality) and lower environmental impact than traditional business models [START_REF] Goedkoop | Product Service-Systems, ecological and economic basics[END_REF][3][4] [START_REF] Wong | Implementation of innovative product service-systems in the consumer goods industry[END_REF] [START_REF] Brandstotter | IT on demand -towards an environmental conscious service system for Vienna (AT)[END_REF]. Product-service paradigm requires the integration of autonomous, geographically distributed and heterogeneous stakeholders in business ecosystem creating, sharing and re-using information across functional teams and enterprise boundaries. In a DPSBE, information management is a fundamental process in order to support after-sales services and next round of product-service innovations to the market. This management is carried out across the boundaries of potentially many enterprises -actors within one DPSBE.

In terms of the collaboration aspects of DPSBEs, they are similar to the so-called Virtual Breeding Environments (VBE) -the pools of organizations and related supporting institutions that have both the potential and the will to cooperate with each other through the establishment of a "base" long-term cooperation agreement and interoperable infrastructure [START_REF] Sánchez | Virtual Breeding Environment: A First Approach to Understanding Working and Sharing Principles[END_REF]. Then, this infrastructure is capable to give birth to so-called virtual enterprises -the temporary networks of independent enterprises (organizations), who join together quickly to exploit fast-changing opportunities and then dissolve [START_REF] Browne | Extended and virtual enterprises -similarities and differences[END_REF].

Research on PSS involves a diversity of scientific and research topics: interoperability of Enterprise Information Systems (EIS), organizational theory, business process management, enterprise architecture, ontology, and knowledge management. Most of these topics are separately largely addressed, but the semantic interoperability of EISs lacks of extended research results in PSS context. Although many EU projects tried to setup a common unified view of business ecosystems, they did not propose any solutions to the problems related to the complex semantic issues underlying interoperability barriers among collaborative enterprises.

The ability of DPSBEs members' information systems to interoperate is one of the key issues in a Product-Service System (PSS). Many projects that deal with interoperability issues, such as ATHENA, INTEROP NoE, ENSEMBLE, UNITE and COIN use semantics as a facilitator for interoperability.

The HORIZON 2020 Framework Programme for research and innovation emphasizes that semantic interoperability is one of the pillars the future research roadmaps. It is important to consider that semantic interoperability goes beyond mere data exchange; it assumes fewer technical preconditions and it deals particularly with data interpretation, where this interpretation is highly dependent on a given context. In this position paper, semantic interoperability is tackled within the perspective of the product-service systems.

So far, none single initiative has addressed the problem of semantic interoperability in PSS context in an emerging Service Economy, to such extent as to manage delivering practical solutions or a larger engagement of the scientific community. Although a range of tools and methodologies for designing PSS exist, they typically represent minor developments of the conventional product design processes, with serious lack of evidence and justification of such tools and methods [START_REF] Baines | State-of-the-art in product service-systems[END_REF].

DISCUSSION

Future research on the semantic interoperability for DPSBEs should build upon the following complementary scientific topics: a) Foundational aspects of semantic interoperability; b) Formal models for DPSBEs; c) Semantically interoperable EISsimplementation and application issues.

The future research on these topics may be considered as a transversal multi-domain approach. Areas to be considered are: complex systems, knowledge management, enterprise interoperability, productservice systems and networked organizations.

Foundational aspects of semantic interoperability

Foundational aspects should deal with theoretical aspects of semantic interoperability -definition, evaluation, knowledge representation requirements and languages.

In defining the foundational aspects for semantic interoperability, we follow the approach of John Sowa [START_REF] Sowa | Knowledge Representation : Logical, Philosophical, and Computational Foundations[END_REF] and define semantic interoperability as follows:

"A sender's system S is semantically operable with a receiver's system R if and only if the following condition holds for any data p that is transmitted from S to R: For every statement q that is implied by p on the system S, there is a statement q' on the system R that: (1) is implied by p on the system R, and (2) is logically equivalent to q. The receiver must at least be able to derive a logically equivalent implication for every implication of the sender's system.".

In semantically interoperable systems, there is no need for any kind of data structures or meta-information which is typically used to assign values so the receiving system can understand the meaning of those values. Instead, exchanged information is considered as a logical statement or a set of the logical statements which describe the semantics of the message from one system to another.

Fig.1. Semantic interoperability of systems

These messages are interpreted by the semantic infrastructure [START_REF] Zdravković | An approach for formalising the supply chain operations[END_REF] (see Fig. 1) which encloses: 1) explicit formal models of the systems (S i ) that need to interoperate (namely, local ontologies O Li ); 2) domain formal models with explicitly defined common semantics (namely, domain ontologies O Di ); and 3) correspondences (mappings) between the concepts of each of the local and domain ontologies (M LiDi ). Then, the mappings between the concepts of two interoperating local ontologies (and corresponding systems' models) can be inferred as logical functions of their mappings with domain ontology.

Formal models for DPSBEs

Formal models for DPSBEs should provide explicit and formal descriptions of the organizational environment whose actors need to interoperate -formal models of enterprise architectures and reference collaboration models. Although the problem of formalizing the enterprise environment is widely addressed by the scientific community, now it has to be reconsidered in context of PSS.

The research of formal models for DPSBEs needs to take into account existing standards, enterprise architectures and models and formal models. Relevant standards include: ISO 15704 -Requirements for Enterprise Reference Architecture and Methodologies [START_REF]ISO 15704, Industrial Automation Systems -Requirements for Enterprise-reference Architectures and Methodologies[END_REF], EN/ISO 19439 -Enterprise Integration -Framework for Enterprise Modeling [START_REF]EN/ISO I9439, Enterprise Integration-Framework for Enterprise Modelling[END_REF] and IEEE 1471 -Recommended Practice for Architectural Description of Software-Intensive Systems-Description [START_REF]IEEE 1471, Recommended Practice for Architectural Description of Software-Intensive Systems[END_REF]. The relevant architectures and models are: Computer Integrated Manufacturing Open System Architecture [START_REF] Amice | CIMOSA -Open System Architecture for CIM[END_REF] (CIMOSA), ARIS [START_REF] Scheer | Business Process Engineering. Reference Models for Industrial Enterprises[END_REF], Zachman Framework [START_REF] Zachman | The Framework for Enterprise Architecture: Background, Description and Utility, Zachman Institute for Advancement[END_REF], The Open Group Architecture Framework [START_REF]TOGAF: The Open Group Architecture Framework[END_REF] (TOGAF), GRAI Integrated Methodology [START_REF] Chen | The GRAI-GIM reference model, architecture and methodology[END_REF],

Purdue Enterprise-Reference Architecture [START_REF] Williams | The Purdue enterprise reference architecture[END_REF] (PERA) and Department of Defense Architecture Framework [START_REF]Definitions and Guidelines[END_REF] (DoDAF). Finally, the formal models for DPSBEs should benefit from the existing work in enterprise ontologies, such as: TOVE (TOronto Virtual Enterprise) ontology [START_REF] Fox | An organization ontology for enterprise modelling: preliminary concepts for linking structure and behaviour[END_REF], The Enterprise Ontology [START_REF] Uschold | The enterprise ontology[END_REF], IDEON TM ontology [START_REF] Madni | IDEONTM: An Extensible Ontology for Designing, Integrating and Managing Collaborative Distributed Enterprises[END_REF] and Supply Chain Operations Reference (SCOR) Ontology [START_REF] Zdravković | An approach for formalising the supply chain operations[END_REF].

Although many of the reference works are complementary, the synthesis of the above results in context of DBPSPEs should take into account a certain amount of existing redundancies and consequently, a need for harmonization. This harmonization has been already started in the development of unifying models, such as Generalized Enterprise-Reference Architecture and Methodology [START_REF]IFAC-IFIP Task Force, GERAM: Generalized Enterprise Reference Architecture and Methodology, Version 1.6.3, IFAC-IFIP Task Force on Architecture for Enterprise Integration[END_REF] (GERAM).

Enterprise architectures developed in the past are contextual, in the sense that they reflect the background and purpose of their developers: CIMOSA for computer integrated manufacturing, GRAI for production management, PERA for system engineering, Zachman for information systems and DoDAF for military operations management. GERAM is considered as the best candidate as a reference architecture to which the concepts of these architectures can be mapped, analyzed and compared [START_REF] Chen | Architectures for enterprise integration and interoperability: Past, present and future[END_REF].

The context in which the above models need to be reconsidered is provided by the trend of "servitization" of products [START_REF] Ducq | Principles of Servitization and Definition of an Architecture for Model Driven Service System Engineering. Enterprise Interoperability[END_REF] and "productization" of services, where these two trends are converging to have a product and a service as a single offering. The occurrence of such PSS have many consequences on the relationships between the manufacturers (including their suppliers and involved services providers) and a customer, related to the business model of a productservice delivery, governance of the process, ownership of the product, etc.

Fig.2. Structure of PSS

The PSS can take a form [START_REF] Tukker | Eight types of Product-Service System: Eight ways to sustainability? Experiences from SUSPRONET[END_REF] (see Fig. 2) of the product-oriented PSS (where, e.g. after-sale services, such as maintenance, repair or recycle are also included in the original act of product sale), use-oriented PSS (e.g. providing the availability of the product, instead of a product itself) or result-oriented PSS (providing the product capability instead of the product).

Obviously, the consideration of the product value chain is different in a PSS than in a conventional productcentric environment, because of the increasing service value element, effect to the customization as well as value of use (including reduced costs due to using less energy or even due to releasing from the responsibilities of asset ownership). This is especially the case in the latter two forms of PSS.

The infrastructure for supporting above scenarios should be provided by the DPSBEs. In a way, DPSBE concept inherits the concept of Virtual Breeding Environment (VBE) for virtual enterprises; it builds upon VBE by considering also new or revised business models, delivery processes, product development processes, etc. It also needs to consider early involvement of the customer in product-service design.

Implementation and application issues of semantically interoperable EISs

Implementation and application issues of semantically interoperable systems are related to the practical aspects -implementation of semantically interoperable EISs, impact to system architectures, reverse engineering methods, etc.

The implementation is based on the semantic infrastructure (see Fig. 3) which includes:

1) formal and explicit domain models (see section 2.2), namely domain ontologies (D oi );

2) formal explicit system models, which transforms the implicit models of the EISs (S i ), used by product manufacturers (PM i ) or service providers (SP i ), to computable formal and explicit ontologies, namely local ontologies (L oi ); and

3) application models, namely problem ontologies (P oi ), which formalize specific problems related to the DPSBEs infrastructure, such as: selection of product or service providers or calculation of the cost of PSS use, rent or lease.

Fig. 3. Semantic infrastructure for DPSBE All three types of the above described formal models in the semantic infrastructure are interrelated by defining the logical correspondences between their respective concepts and relations, so to achieve a semantic interoperability of the systems which are using those (see section 2.1). These mappings are stored in unifying or mapping ontology (M o ).

Then, access to the overall knowledge of DPSBE (knowledge-of-systems and knowledge-of-domain) is enabled through the collaborative semantic applications (SA i ), used by the coordination service providers of DPSBE (CP i ).

Each of the semantic application is used to resolve individual problem of DPSBE, formalized by the problem ontologies (P oi ) which import the central mapping ontology.

CONCLUSIONS

The future work on pursuing the topics above should be strongly inter-related, because it involves the research on the three different perspectives on the individual problem: modelling, organizational and technical. At the same time, it addresses theoretical, foundational and practical issues. Thus, it is expected to provide the holistic view on the semantic interoperability in DPSBEs.

In this paper, we propose the foundations of semantic interoperability, candidate formal models for DPSBEs and approach to develop a semantic framework which would achieve the semantic interoperability.

The high innovative potential of the results of the research above draws from the fact that it enables organizations to interoperate without a need to fulfill the complex technical requirements. Thus, interoperation may easily become accessible (even as a service) for SMEs or any other organization that wants to participate in DPSBEs. Obviously, capability to quickly establish an infrastructure supporting the specific PSS is critical, since it contributes to the increased openness and flexibility and decreased interaction costs in DPSBE.

The end-users of the DPSBEs research are all organisations engaging in collaborative activities, especially PSS actions, whose performance depends on their capability to exchange information to support after-sales services and next round of product-service innovations to the market. One of the success factors of the industrial SMEs is their ability to participate in as many supply chains as possible, where each of them poses different manufacturing or service requirements and standards. Semantic interoperability could help in transferring this diversity manageable. Therefore, it is expected to increase the flexibility of SMEs and consequently, their capability to take part in DPSBEs.

In the long term, the advances in research on DPSBEs could lead to the establishment of new PSSs, in a global market undergoing a transition from being a product provider into being a value provider by integrating services into core-products offerings and vice-versa.

ACKNOWLEDGEMENT

The work presented in this paper was supported by the program for scientific cooperation between Serbia and France PHC PAVLE SAVIC (project No. 23494VF 2010/2011); and the Ministry of Education and Science of the Republic of Serbia (project No. III41017).


Guillaume Melquiond
email: guillaume.melquiond@inria.fr

Floating-point arithmetic in the Coq system

Keywords: oating-point arithmetic, formal proofs, Coq system

The process of proving some mathematical theorems can be greatly reduced by relying on numerically-intensive computations with a certied arithmetic. This article presents a formalization of oating-point arithmetic that makes it possible to eciently compute inside the proofs of the Coq system. This certied library is a multi-radix and multi-precision implementation free from underow and overow. It provides the basic arithmetic operators and a few elementary functions.

Introduction

Some mathematical theorems have recently been proved by performing huge amounts of numerical computations (e.g. Hales' proof of Kepler's conjecture), hence casting some doubts on the validity of their proofs. By performing these numerical computations inside a formal system, a much higher condence in the theorems would have been achieved. Indeed, these systems ruthlessly track unproven assumptions and incorrectly applied theorems, hence helping the user to perform foolproof reasoning. Fortunately, in some of these systems, especially in those based on type-theoretical formalisms, the ability to perform calculations is steadily increasing. It then becomes possible to use these calculations in place of some traditional deductive proofs, hence freeing the user from this burden [START_REF] Gonthier | Formal proof the four-color theorem[END_REF].

One of these formal systems is the Coq proof assistant [START_REF] Bertot | Interactive Theorem Proving and Program Development. Coq'Art: the Calculus of Inductive Constructions[END_REF], which is based on the calculus of inductive constructions. Its formalism makes it possible to evaluate functions and to use their results inside proofs. This system is therefore a good candidate for implementing certied yet ecient arithmetics, and hence for using numerical computations to formally prove mathematical results [START_REF] Grégoire | A purely functional library for modular arithmetic and its application to certifying large prime numbers[END_REF].

Irrespective of these considerations, oating-point arithmetic is widely used in computer programs as a fast approximation of the traditional arithmetic on real numbers. By design (limited range, limited precision), it is ecient for performing numerically-intensive calculations. While these calculations are most often encountered in numerical simulations of physical phenomenons, Thomas Hales has shown that they could also be a great help for proving some mathematical theorems [START_REF] Hales | A proof of the Kepler conjecture[END_REF].

This paper presents a formalization of an eective oating-point arithmetic inside the Coq system. The incentive is to provide a library that can be used to prove theorems by performing numerical computations. For instance, interval arithmetic is a simple method for computing enclosures of real-valued expressions and therefore to prove bounds on such expressions. Moreover, verifying such bounds is usually no harder than computing them directly, so a formal system would end up having to do these computations, even when helped by an external oracle. The oating-point library presented in this paper provides the denitions and the corresponding correctness lemmas that make it easy to write tools that performs proofs by numerical computations [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF].

A formal system is an unusual environment for developing a oating-point arithmetic, so Section 2 presents some design decisions and how this library compares to other oating-point formalizations. Section 3 then details the implementation of the basic operators: addition, multiplication, division, and square root. The library also encompasses some elementary functions which are described in Section 4. Finally, some recent realizations and future works are presented in Section 5.

Floating-point arithmetic

Floating-point numbers are usually a subset of rational numbers, with some additional values for handling exceptions (e.g. innities). A radix β is associated to a oating-point arithmetic, and its nite numbers can be represented as m•β e , with m and e two integers. Most common radixes are β = 2, widely available in general-purpose processors, and β = 10, often found in nancial applications. One can also nd β = 16 in some older hardware, and β = 3 in exotic computers.

For the sake of realizability and/or space eciency, a precision is usually set, that is an integer p such that the oating-point numbers are restricted to mantissas m that are bounded: |m| < β p . For the same reasons, exponents e are also constrained to a range [e min , e max ]. For instance, the double precision arithmetic described in the IEEE-754 standard [START_REF] Stevenson | An American national standard: IEEE standard for binary oating point arithmetic[END_REF] is a radix-2 arithmetic with p = 53, e min = -1074, and e max = 971. A multi-precision library like MPFR [START_REF] Fousse | MPFR: A multiple-precision binary oating-point library with correct rounding[END_REF] works with any precision but still has bounded exponents, though the bounds are suciently big so that they do not matter usually.

Number format

The oating-point formalization presented in this article supports any radix β ≥ 2. Indeed, because of their simplicity, the operations described in Section 3 do not rely on the properties of a specic radix. While any radix can be used, the library does not handle mixed-radix operations. One cannot add a radix-2 number with a radix-3 number and obtain a radix-10 number. So the radix can be seen as a global setting for the oating-point arithmetic. This is not the case for the precision of the numbers. It can be set independently for each operation, and an operation will return a number with the given precision.

Be it in Coq or in other systems, integers are mathematical objects with no immediate notion of limited range. In particular, the integers m and e representing a number m • β e can be arbitrarily big. So there is no need for an upper bound on the precision available to computations. Similarly, constraining the exponents to a bounded range would be articial. So the formalization has neither e min nor e max . As a consequence, underow and overow no longer occur during computations. It means that exceptional values like signed zeros, subnormal numbers, and innities, are a lot less useful, and have hence been discarded from the formalization.

The unbounded range of exponents has some immediate properties. First, a oating-point operation will only return zero when the ideal mathematical result is zero too. Second, all the results are normal numbers, so bounds on the relative errors are always veried, which makes it easier to write oating-point algorithms and to perform their error analysis. One no longer has to deal with the traditional sentence: The result is correct assuming that no underow nor overow occurred in the course of the computations.

The formalization nonetheless supports an exceptional value: Not-a-Number.

It is returned by oating-point operators, when the corresponding mathematical function on real numbers is not dened on the inputs. For instance, a NaN will be returned when computing 1 0 and √ -1. As usual, this exceptional value is an absorbing element for all the oating-point operations. So the nal result of a sequence of computations will be NaN, if any of the intermediate computations had invalid inputs. This is especially useful, since the pure functional programing language of Coq does not oer native exceptions and traps, so delaying the handling of the invalid cases simplies the use of oating-point operators inside proofs.

Data sets and functions

Since neither the precision nor the exponent are bounded, any number m•β e can be represented and handled within this formalization. Let us note F β this subset { m • β e | (m, e) ∈ Z 2 } of the real numbers. These numbers will be represented as pairs of integers (m, e). Notice that these pairs are not normalized a priori, so (m, e) and (m • β, e -1) are two dierent representations of the same real number. The set F β will denote the whole set of oating-point numbers, that is F β extended with a NaN value.

To summarize, a oating-point operation like the division will be a function with a signature depending on its (implicit) rst parameter β:

Fdiv : ∀β : radix, rounding → precision → F β → F β → F β
The rounding type contains modes for selecting a oating-point result when the ideal mathematical value is either outside of F β or not representable with the required precision. The supported modes are detailed in Section 3.1. The precision type denotes all the positive integers, though the results will not be specied when precision is 1.

So that the oating-point division can be used in a meaningful way, the library contains the following Coq theorem: For any two oating-point numbers

x and y in radix β (bigger than 1) and for any rounding mode and precision prec, the result of the algorithm represents the same real number as the real quotient

x y once rounded to a oating-point number.

Theorem Fdiv_correct : forall radix mode prec ( x y : float radix ) , 1 < radix -> Fdiv mode prec x y = round radix mode prec ( x / y ).

This theorem matches the IEEE-754 requirement: Every operation shall be performed as if it rst produced an intermediate result correct to innite precision and with unbounded range, and then rounded that result according to. . . Note that the real quotient x y and its rounded value are dened in a purely axiomatic way in Coq: Standard real numbers are an abstract type, so they are unusable for computing. All the computations are performed on integers and they are proved afterwards to match the abstract result.

Computability

Several formalizations of oating-point arithmetic exist, for various proof assistants, and they have been successful in proving numerous facts on oatingpoint algorithms. For the sake of conciseness, only three of them will be cited here [START_REF] Russino | A mechanically checked proof of IEEE compliance of the AMD-K7 oating point multiplication, division, and square root instructions[END_REF][START_REF] Harrison | A machine-checked theory of oating-point arithmetic[END_REF][START_REF] Daumas | A generic library of oating-point numbers and its application to exact computing[END_REF]. The main dierence between these libraries and the one described in this paper is in the rationale. They were designed for proving oating-point code, but not for performing actual computations. For instance, the following predicate is the usual characterization of rounding to -∞ in their formalization: A real number x is rounded to f in a oating-point format (radix, precision, and exponent range) if

f ∈ format ∧ f ≤ x ∧ ∀g ∈ format, g ≤ x ⇒ g ≤ f.
While this is the best denition for describing the behavior of a oatingpoint rounding, this predicate does not provide any computational content, since there is no direct way for computing f given this denition and x. For their verication purpose, these libraries do not need to provide explicit algorithms for rounding values. If they do (either directly or through an oracle), the motivation is to evaluate oating-point operations only for particular explicit values [START_REF] Harrison | A machine-checked theory of oating-point arithmetic[END_REF]. They were not designed to perform numerically-intensive computations involving millions of oating-point operations, as is done in the example of Section 5.1 [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF].

Note that the theorems provided by the existing Coq libraries [START_REF] Daumas | A generic library of oating-point numbers and its application to exact computing[END_REF] cannot be easily reused, unfortunately. Indeed, since their goal is the verication of code relying on actual oating-point arithmetic, e.g. IEEE-754, the bounded exponent range (especially the minimal exponent) is intimately embedded in their formalization.

Modules

The rst part of the formalization denes and proves algorithms that work for any radix. While usable to perform oating-point computations, the purpose of this generic implementation is more of a reference implementation. For instance, counting the number of digits of a mantissa m is performed by increasing n until β n is bigger than |m|. This is correct yet inecient. The issue is similar when shifting mantissas left m•β n and right ⌊m/β n ⌋ since the shifts are performed by powering followed by a multiplication or a division.

The Z integers provided by Coq's standard library are represented as list of bits. As a consequence, when the radix is β = 2, mantissa operations can be performed much more eciently. For instance, shifts can be computed by adding zeros at the head of the list or removing the rst elements, since the least-signicant bits are stored rst. So the library has been specialized for specic radixes in order to improve performances.

This specialization is achieved thanks to the module system of Coq. First of all, an interface (module type ) FloatOps contains the signature and the specication of all the oating-point operations. A module that implements this interface must therefore provide a set of oating-point operators and a proof that these functions return the exact same results as the functions of the reference implementation. For instance, if the module M has type FloatOps, the M.add function is a oating-point adder and the M.add_correct theorem proves its correctness. An important point is that the radix is no longer a parameter of the oating-point operations.

The rst implementation of this interface is the module functor GenericOps. Given a radix, it generates a module of type FloatOps which is a simple wrapper around the reference implementation. For instance, Module M := GenericOps Radix10 generates a module M which provides decimal oating-point arithmetic. The second implementation is provided by the functor SpecificOps. Given a set of functions for performing special operations on mantissas (e.g. shifts), this functor generates an improved module that directly manipulates integers at the digit level. For instance, SpecificOps StdZRadix2 takes advantage of the representation of Coq's standard integers as bit lists in order to speed up radix-2 oating-point computations.

The implementation can be sped up even further by replacing the bit lists with binary trees of xed-size integers [START_REF] Grégoire | A purely functional library for modular arithmetic and its application to certifying large prime numbers[END_REF]. The binary-tree structure allows for divide-and-conquer algorithms, e.g. Karatsuba's multiplication. Moreover, Coq can use 31-bit machine integers for representing the leaves of the trees. This considerably reduces the memory footprint of the integers, and makes it possible to delegate arithmetic computations on the leaves to the processor. Arithmetic operations are no longer performed one bit after the other, but by blocks of 31 bits at once. As a consequence, the module SpecificOps BigIntRadix2 is currently the fastest radix-2 implementation, when evaluating expressions with Coq's virtual machine.

Basic operators

This section presents the implementation of the basic arithmetic operators.

These operators have been developed with a simple design on purpose, so that they can be easily proved correct. As a consequence, this Coq library does not follow the laudable philosophy of the multi-precision oating-point library MPFR: The complexity should, while preserving correct rounding, depend on the precision required on the result rather than on the precision given on the operands [START_REF] Fousse | MPFR: A multiple-precision binary oating-point library with correct rounding[END_REF]. When working at constant precision, this issue fortunately does not matter.

Rounding

If the result x of an exact operation can be represented as a pair (m, e) with |m| < β p , then this oating-point value should be returned by the associated operator running at precision p. Otherwise, a oating-point number near x is chosen according to a rounding direction and is returned by the operator.

Let us assume that x is a positive real number. Let e = ⌊log β x⌋ -p + 1 and m = ⌊x • β -e ⌋. Both m • β e and (m + 1) • β e can be represented by oating-point numbers with mantissas of p digits. Moreover, they are respectively the biggest number smaller than x and the smallest number bigger than x. These are the two candidates toward which x can be rounded.

We rst need to have the position of x relatively to these two numbers. Let us pose d = x • β -e -m. By construction, d is in the range [0, 1). If d is 0, then x is equal to m • β e . This position is called pos_Eq in the library. Otherwise, d is either smaller than 1 2 (pos_Lo), equal to 1 2 (pos_Mi, x is half-way between m • β e and (m + 1) • β e ), or bigger (pos_Hi).

Choosing the correctly-rounded value is therefore a simple case study. For instance, when rounding to nearest the number x at position pos_Hi, the oatingpoint value (m+1, e) is returned. This four-position system is strictly equivalent to the usual hardware approach of using two bits: the rounding bit and the sticky bit [START_REF] Ercegovac | Digital Arithmetic[END_REF]. The four rounding directions of the IEEE-754 standard are supported: toward zero, toward -∞, toward +∞, and to nearest (with tie breaking to numbers with even mantissa). New rounding modes can easily be supported, as long the discontinuity points are either oating-point numbers or half-way between consecutive ones.

The position system does not have to be restricted to powers of the radix though. Its generalization will give the correctness proof of the division operator for free (Section 3.3). The library denes a correctly_located predicate to express that the mantissa m and the position pos are correctly chosen for a given scale. This predicate holds when:

   x = m × scale if pos is pos_Eq m × scale < x < (m + 1 2 ) × scale if pos is pos_Lo . . .
For arbitrary numbers x and scale, the position and m are not computable. 2But they are when we already know an integer m ′ and a position pos ′ for a scale scale/n, with n a positive integer. Indeed, m is then ⌊m ′ /n⌋, and the position is given by the adjust_pos function of the library, which compares the remainder m ′ -n × ⌊m ′ /n⌋ with ⌊ n 2 ⌋. For instance, the new position is pos_Mi if the remainder is n-1 2 and the old position is pos_Mi too.

Addition and multiplication

The set F β is actually a ring for the addition and the multiplication inherited from the real numbers. For instance, given two numbers (m 1 , e 1 ) and (m 2 , e 2 ), the pair (m 1 • m 2 , e 1 + e 2 ) represents their product. For the addition, the mantissas have rst to be aligned so that the two numbers are represented with the same exponent. The mantissas can then be added. Functions for exact addition and exact multiplication are therefore provided. This ring structure also helps in dening and proving the rounded operators.

Indeed, if the mantissa of the exact result (m, e) has less than p digits, this is also the correctly-rounded result. Otherwise, (m, pos_Eq) correctly locates the number m • β e with the scale β e . So it can be scaled down to the lower precision p by using the algorithm of Section 3.1. Once done, the new position is used to decide which number to return.

Notice that, for the oating-point rounded addition, this approach is especially inecient when the exponents e 1 and e 2 are far from each other, as it requires a huge shift followed by a long addition. For instance, the sum of 1 • β 10000 and 1 • β 0 involves a 10001-digit addition, while the result is more or less trivial for small precisions, depending on the rounding direction. A better approach would be to extract from the inputs the mantissa parts that actually matter in the rounded result.

Division and square root

For division and square root, one cannot expect to compute the exact values rst. It is nonetheless possible to perform divisions and square roots on integer mantissas and to obtain an exact integer remainder.

Given two positive numbers (m 1 , e 1 ) and (m 2 , e 2 ), the division operator rst computes the lengths of the mantissas l 1 = 1+⌊log β m 1 ⌋ and l 2 = 1+⌊log β m 2 ⌋. The integer n = max(l 2 + p -l 1 , 0) is such that the integer quotient q = ⌊m 1 • β n /m 2 ⌋ has at least p digits. The exact result m 1 /m 2 • β e1-e2 of the division is correctly located by (m 1 • β n , pos_Eq) with the scale β e1-e2-n /m 2 . So it can rst be scaled down by a factor m 2 , which gives a new location with q as a mantissa. It can then be scaled again by a factor β k if q has p + k digits. In the end, we get a location of x whose scale is a power of β and whose mantissa has p digits exactly. This is sucient to get the correctly-rounded value of x.

The square root algorithm cannot rely on these changes of scales. It is nonetheless similar. Indeed, in Section 3.1, the rounding algorithm relies on the remainder of the euclidean division. So it can be adapted so that it uses the remainder of the integer square root. In order to compute the rounded result of the square root of (m, e), the operator rst computes the length l of m. It then chooses the rst integer n bigger than max(2p -l -1, 0) and such that e -n is an even integer. The integer square root s = ⌊ √ m • β n ⌋ has at least p digits. The remainder is r = m • β n -s 2 . If r is zero, the exact result is at the (s, pos_Eq) location with a scale β (e-n)/2 . Otherwise, the position is obtained by comparing r to s + 1 2 , which is the half-distance between s 2 and (s + 1) 2 . Since r and s are both integer, they can actually be compared directly. Finally, the location is scaled down again so that the mantissa has exactly p digits.

Elementary functions

For the basic operations, the exact real result can either be represented directly as a oating-point number, or with the help of a representable remainder. This is no longer the case for elementary functions. Except for a few trivial inputs, e.g. 0, one can only compute non-singleton ranges enclosing the exact result. This is nonetheless sucient in order to get correct rounding, as shown by Ziv's iterative process [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF].

Formalizing this process in Coq, however, depends on theorems that are currently out of scope. So the elementary functions do not return the correctlyrounded result. Instead, they return an interval enclosing the exact mathematical result. Fortunately, this interval is sucient when proving inequalities on real-valued expressions. So the Coq tactics [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF] that depend on this oating-point library do not need the correct rounding.

The library currently supports the functions cos, sin, tan, arctan, and exp. The implementation of other elementary functions like log is a work in progress. The library relies on an exact division by 2 for oating-point numbers, so the elementary functions cannot be used when the radix is odd.

Series evaluation and interval arithmetic

Elementary functions are evaluated thanks to simple power series whose terms happen to be alternating and decreasing for small inputs (say 0 ≤ x 2 ≤ 1 4 ). For instance,

arctan x = x • ∞ i=0 (-1) i • (x 2 ) i 2i + 1
From |x| ≤ 1 2 , we can decide how many terms are needed so that the remainder of the series is guaranteed to be smaller than a threshold β -k . For instance, when β ≤ 4, k terms are sucient. This does not mean that the function will compute that many terms, it just means that k is the explicitly-constructed decreasing argument needed to dene the recursive summation. Actually, the function tests the current term against the threshold and stops as soon as it is smaller.

A careful error analysis would then permit us to dene and to prove an algorithm for evaluating truncated summation with an absolute error less than β -k too. The relative error for arctan would be similar when computing it this way. Unfortunately, this error analysis is currently out of scope. As a replacement, a small kernel of oating-point interval arithmetic was implemented and proved in Coq. So the summation of the series is instead performed with intervals at precision k. This takes into account both the rounding errors and the series remainder, and it trivially ensures that arctan x is contained in the computed range. But the relative width of the range is no longer guaranteed to be smaller than β -k and hence to converge toward zero when k increases. This prevents the completion of a proof that Ziv's process actually terminates for this implementation. In order to get an input x with an absolute value smaller than 1 2 , an argument reduction is performed. For the three direct trigonometric functions, anglehalving formulas are used:

cos(2 • x) = 2 • (cos x) 2 -1 sign(sin(2 • x)) = sign(sin x) • sign(cos x)
These formulas give cos x and the sign of sin x for any x. The values sin x and tan x are then reconstructed thanks to the following formulas:

sin x = sign(sin x) • 1 -(cos x) 2 tan x = sign(sin x) • sign(cos x) • (cos x) -2 -1
When no argument reduction is needed for sin and tan, the library does not rely on the power series of cos. Instead, it relies directly on the series of sin.

For the exp function, the argument is rst reduced to a negative number by using exp x = (exp(-x)) -1 if needed. The argument is then brought in the domain [- 1 2 , 0] by squaring the result: exp x = (exp x 2 ) 2 . On this reduced domain, the power series of exp is alternated, hence computable the same way than the other functions.

These argument reductions were chosen because they provide exact reduced arguments and are simple to prove. It would be interesting to compare their performances against the more traditional approach of using a trivial reconstruction and inexact reduced arguments:

cos x = cos(x -k • 2π) with k = x 2π exp x = exp(x -k • log β) • β k with k = x log β 4.2.

arctan and π

The argument of the arctan function is rst reduced to x > 0 by using the parity of the function. The result is then computed thanks to the power series on the domain [-1 3 , 1 2 ], after a potential reduction with the following formulas:

arctan x = π 4 + arctan x-1 x+1 for x ∈ [ 1 2 , 2] π 2 -arctan 1 x for x ≥ 2
Notice that the result reconstruction of arctan involves the constant π 4 . It is computed thanks to Machin's formula 4 • arctan 1 5 -arctan 1 239 , which needs the computation of arctan for small arguments only. In order not to recompute this constant for every arctan evaluation, it is stored inside a co-inductive object. The i-th element of the co-inductive object is dened as an interval computed at precision 31 • i with Machin's formula and hence containing π 4 . A co-inductive object can be seen as an innite list. One can destroy it into two elements: its head and its tail which is again an innite list. Whenever a co-inductive object is destructed into these two elements, Coq's virtual machine has to compute the value of its head. 3 Moreover, Coq's virtual machine remembers this head value, so that it can be instantly reused the next time a function destructs the co-inductive object. Therefore, co-inductive objects can be used as a cache (that would never be ushed though). In particular, the interval constant for π at a given precision (with a granularity of 31 digits) is computed only once per Coq session.

Reconstruction and accuracy

The library performs all the intermediate computations at the same precision. The reconstruction process may, however, incur a loss of accuracy. For instance, in the case of exp, one bit (or log β 2 digits) of accuracy is lost for each squaring, as long as second-order rounding errors are negligible. Figure 1 shows the number of bits lost (relative width of the output interval on a log 2 scale) depending on the magnitude of the input (again on a log 2 scale).

For each input exponent, the relative width was averaged on 128 values. The precision used for the computations was 30 bits, but the overall gure does not change with other precisions: Be it 50 or 100 bits, the plot looks similar. Figure 1 shows that the loss is constant as long as the input is less than 2 -1 . Then one additional bit is lost each time the input gets twice as big, except for the arctan function, whose reconstruction process has a xed number of operations.

As a consequence, the functions have been modied to take into account the loss of accuracy caused by the argument reduction: The internal precision is increased depending on the magnitude of the input, so that the resulting interval has an average width of one unit in the last place. Therefore, even if the elementary functions do not provide correct rounding, they still have an accuracy similar to the basic arithmetic operators.

Time complexity

Figure 2 shows the time it takes to Coq for computing the result of an elementary function (β = 2) depending on the internal precision. Note that the input value is reduced to 0.4 which will be used while evaluating the power series. The cos plot is identical to the sin one: The cost of the square root while rebuilding sin is negligible with respect to the other operations.

The last three plots are for the arctan function. Indeed, since the input value is big, the computations will have to access the cache containing enclosures of stored: binary trees with 31-bit native integers at the leaves. Indeed, the steps occur around precisions 31 • 2 5 = 992, 31 • 2 6 = 1984, 31 • 2 7 , and so on.

The slope of the exp and sin plots is about 2.4, which means that the time complexity of these functions is asymptotically p 2.4 .

MPFR test data

The elementary functions of this library have been compared to the corresponding MPFR radix-2 functions. The MPFR library is shipped with test les containing input values and the expected results for these functions at a precision of 53 bits. All these input values were tested with the Coq implementation in order to compare it with MPFR performance-wise.

The strategy for computing the results in Coq is as follows. An interval enclosing the exact result is computed with an internal precision p i . If both bounds round to the same oating-point number at precision 53, then this number is the correctly-rounded result. Otherwise the interval is too wide and the computation starts again at precision p i+1 = ⌊p i × 3 2 ⌋. And so on until a result is returned. The starting precision p 1 = 63 is arbitrarily chosen, as is the 3 2 multiplier. Note that this Ziv strategy does not require any termination proof: Since the input set is xed, a suitable upper bound on i can be found by experimenting.

Table 1 gives the average slowdown caused by using Coq instead of MPFR to compute correctly-rounded results for the test data. Note that these results are only meant to give a rough idea of the relative performance of the implementations, since the evaluation strategy is arbitrary and the input values are unusually large. The table shows the slowdown depending on the integer type the library is instantiated with. The rst column is obtained for fast binary-tree integers, while the second column is obtained for standard bit-list integers. First noticeable point: Computing a value with Coq is about 2000 times slower than computing it with an optimized C program. This is the expected slowdown magnitude when programming inside the formal system. It is hopefully acceptable for proving some theorems with numerical computations.

Second point: Using standard integers only incurs an additional slowdown of 7× for arctan and cos with respect to fast integers. This one is a bit puzzling: Since integer computations are much faster on the processor, the slowdown should be bigger. This may be explained by the small size of the integers. Moreover, some common operations are trivial on standard integers, e.g. shifts.

Third point: There is yet another slowdown when Coq computes sin and tan with standard integers. It is due to the use of a square root when reconstructing the nal result. The square root on standard integers does not only return the integer result, but it also returns a proof that this integer is the correct result. In the type-theoretical setting of Coq, this proof is an explicit λ-term whose body and type depend on the integer input. Coq wastes a long time (65% of the total time) building the whole proof term, while it is irrelevant from a computational point of view. 4 This explains why the binary-tree integers are much faster there: The correctness of their square root is guaranteed by a separate theorem instead of being carried by the code.

These tests also helped to detect an unexpectedly unproven (and wrong!) theorem in the Coq library that formalizes native 31-bit integer computations. 5They also showed that the test values of MPFR were not as dicult as they could be, since the Coq functions could almost always get the correct result at the second iteration, hence with an internal precision of 94 bits.

Reducing arguments further

The functions of the library only perform an argument reduction until the input is smaller than 1 2 . This is unusual [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF], since a smaller argument will speed up the power series evaluation. Figure 3 shows the time variation depending on the threshold chosen for the reduction of exp. To account for the longer reconstruction, the precision of the series evaluation was increased by one whenever the threshold was divided by two. The plot shows that reducing an argument until between 2 -7 and 2 -13 is the best choice for MPFR data. When reducing the argument even further, the cost of the reconstruction steps becomes too large, hence negating the speedup of the power series. Therefore, changing the threshold brings the evaluation of exp in the same speed range as the trigonometric functions, about ×2000 slower than MPFR.

Note that this optimal threshold cannot be generalized to other precisions, since the number of operations during the evaluation of the power series depends on the precision, while the reconstruction depends only on the magnitude of the input.

For the trigonometric functions, the speedup when reducing the threshold is less noticeable. This may be due to the small amount of information available in an approximation of cos x for x small. So one could hope to achieve some speedup by evaluating instead a function umc x = 1 -cos x and using the argument reduction

umc(2 • x) = 2 • (umc x) • (2 -umc x).
Unfortunately, it does not.

Conclusion

Proving theorems

Implementing oating-point arithmetic inside a formal system is hardly useful on its own. This arithmetic is only a means to eciently prove mathematical theorems. So some Coq tactics have been developed. They rely on oating-point computations in order to automatically prove bounds on real-valued expressions [START_REF] Melquiond | Proving bounds on real-valued functions with computations[END_REF]. They perform interval arithmetic and their usage is similar to those of some existing PVS strategies [START_REF] Muñoz | Real number calculations and theorem proving[END_REF]. In addition, these tactics support bisection search, and they rely on a bit of automatic dierentiation for doing interval evaluation with Taylor's order-1 decomposition. As a consequence, these Coq tactics are able to automatically handle a theorem originally proved in PVS [START_REF] Daumas | Guaranteed proofs using interval arithmetic[END_REF]. The exact same formal methods are used for proving this theorem, but the Coq tactics do not need to rely on an external oracle.

The proved theorem states a tight bound on the relative error between Earth local radius

r p (φ) = a 1 + (1 -f ) 2 × tan 2 φ
and a degree-10 polynomial with single-precision oating-point coecients that is approximating it. The original PVS proof was composed of about 10000 generated scripts. It took several hours to check all of them on a 48-core machine in 2005. The Coq proof, a 25-line script, took a few minutes to check on a single core in 2008.

While the Coq tactics are performing too many computations because there is no oracle, they benet from the use of oating-point arithmetic. Indeed, PVS' strategies are using rational numbers as interval bounds. As a consequence, the numbers that appear in the intermediate computations of the PVS proof carry thousands of digits, since all the bounds are computed exactly. 6 On the contrary, all the oating-point bounds in the Coq proof are rounded outwards to a precision of 30 bits. So the computations are not slowed down by the size of the numbers, which explains the tremendous speedup that is achieved on this example.

Coq was also used to check the correctness of a single-precision polynomial approximation of exp [START_REF] Tang | Table-driven implementation of the logarithm function in IEEE oating-point arithmetic[END_REF]. Again, it only took a few minutes for the formal system to automatically prove that the approximation is accurate enough.

Future works

The oating-point formalization described in this article is part of the Coq library available at http://www.lri.fr/~melquion/soft/coq-interval/ The long-term goal is to develop and adapt numerical methods for a proof assistant, so that computational proofs like Hales' one can be completely checked inside a formal system. Implementing a oating-point arithmetic inside Coq is a step toward this goal, and there is still space for improvements.

Obviously, new elementary functions should be added, so that the usual ones at least are available. Fortunately, with a suitable argument reduction, functions like log can also be expressed as alternating series. So their implementation and formal proof should closely match the ones for the existing functions, hence making them straightforward.

More importantly, the interval evaluation of elementary functions should be replaced by a static error analysis. There are currently no formal methods for doing this kind of analysis for multi-precision algorithms, so this will rst require to build a comprehensive formalism. Not counting the ability to actually certify multi-precision algorithms, there are two benets to a formalized static analysis. First, removing intervals will speed up the functions a bit. Second, it will allow us to implement Ziv's strategy [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF] and get correctly-rounded results without relying too much on axioms.

While correct rounding is not needed for proofs, it would be a great help in writing MPFR-based oracles. Indeed, since correct rounding allows for portable results, a numerical computation that succeeds in the oracle would also succeed in Coq. As a consequence, it would become possible to carefully craft the oracle so that the proofs it generates need as few computations as possible to prove a given theorem in Coq.

The other way around, having portable results makes it possible to directly perform the extraction of a high-level numerical algorithm written in Coq to a compilable language (e.g. Ocaml) with bindings to the arithmetic operators and elementary functions of MPFR. That way, both the development and the certication of a numerical application could be done in Coq, while its execution would be as fast as currently possible for a multi-precision code.

4. 2 .

 2 Argument reduction 4.2.1. Forward trigonometric functions and exponential

Figure 1 :

 1 Figure 1: Relative width of the result depending on the normalized exponent of the input.

Figure 2 :

 2 Figure 2: Time for computing a function at point 111111111111 • 2 -30 ≃ 103.5.

Figure 3 :

 3 Figure 3: Time needed for computing exp on MPFR data depending on the argumentreduction threshold.

Table 1 :

 1 Slowdown of Coq functions with respect to MPFR functions.

	Function Coq fast ints Coq std ints
	cos	1700×	12000×
	sin	2000×	40000×
	tan	1900×	37000×
	exp	3300×	21000×
	arctan	470×	2000×

More generally, rounding operators are functions that are not continuous, so they are not computable on the set of real numbers.

Since the co-inductive object is an innite list, its elements cannot be computed at the time the object is created; they are lazily evaluated when the object is destructed.

. The plots correspond to various cache behaviors. The rst plot happens with a hot cache: π 4 is already available. The second plot happens with a cold cache: Coq has to ll it. As shown on the gure, initializing the cache takes a huge time at higher precisions. Indeed, since the cache is stored inside a coinductive list, accessing the value at a given precision of the constant requires to initialize all the values at lower precisions. This issue can be avoided by adding an indirection: Instead of directly storing the constants in the list, one can store trivial co-inductive objects that do compute the constants. This method gives the third plot. It is the one nally implemented.Some regular steps are visible on the plots, especially for arctan at higher precisions. These steps are most certainly related to the way fast integers are

By using a lazy evaluation scheme, the proof term would never be computed, since it is not used. It would, however, slow down all the other computations since none of them would benet from lazy computations.

The correctness proof for the

62-bit square root algorithm was not proved in the standard library of Coq, only admitted. In fact, it would have been impossible to prove it, as the Newton iteration used in the algorithm was actually incorrect. This was detected thanks to discrepancies with MPFR results for large arguments to the sin function.

Note that the PVS computations are only exact for truncated power series, not for the complete elementary functions.