Laurent Thévenoux
email: laurent.thevenoux@inria.fr.

Philippe Langlois

Matthieu Martel

Automatic source-to-source error compensation of floating-point programs: code synthesis to optimize accuracy and time

Keywords: code synthesis, compensation, execution-time performance, floating-point arithmetic, multi-criteria optimization, numerical accuracy

Numerical programs with IEEE 754 floating-point computations may suffer from inaccuracies, since finite precision arithmetic is an approximation of real arithmetic. Solutions that reduce the loss of accuracy are available, such as compensated algorithms or double-double precision floating-point arithmetic. Our goal is to automatically improve the numerical quality of a numerical program with the smallest impact on its performance. We define and implement source code transformations in order to derive automatically compensated programs. We present several experimental results to compare the transformed programs and existing solutions. The transformed programs are as accurate and efficient as the implementations of compensated algorithms when the latter exist. Furthermore, we propose some transformation strategies allowing us to partially improve the accuracy of programs and to tune the impact on execution time. Tradeoffs between accuracy and performance are assured by code synthesis. Experimental results show that userdefined trade-offs are achievable in a reasonable amount of time, with the help of the tools we present in the paper.

INTRODUCTION

In this paper, we focus on numerical programs using IEEE 754 floating-point arithmetic. Several techniques have been introduced to improve the accuracy of numerical algorithms, such as expansions [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF][START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF], compensations [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF], differentiation methods [START_REF] Langlois | Automatic linear correction of rounding errors[END_REF] or extended precision arithmetic using multiple-precision libraries [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF][START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF]. Nevertheless, there are numerous and clearly identified numerical failures [START_REF] Blair | Patriot missile defense: Software problem led to system failure at Dhahran, Saudi Arabia[END_REF][START_REF] Lions | Ariane 5 flight 501 failure, report by the inquiry board[END_REF]. This illustrates that these improvement techniques are perhaps not so widely used when writing numerical software, or not sufficiently automated to be applied more systematically. For example, the programmer has to modify the source code by overloading floatingpoint types with double-double arithmetic [START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF] or, less easily, by compensating the floating-point operations with error-free transformations (EFT) [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF]. The latter transformations are difficult to implement without the preliminary step of designing the modified algorithm, which is typically done manually.

As much as the numerical accuracy, execution time is a critical matter for programs, especially in embedded systems where reducing the execution time may lead to less energy consumption or better reactivity. Taking into account both accuracy and execution time simultaneously is a difficult problem because these two criteria do not cohabit well: improving accuracy may be costly and execution-time improvement can impact accuracy as well. A solution is to provide trade-offs between accuracy and time [START_REF] Langlois | Automatic code transformation to optimize accuracy and speed in floatingpoint arithmetic[END_REF].

We present a method that allows a usual developer to automatically improve the numerical accuracy of his programs without increasing execution times significantly. Although we do not provide error bounds on the processed algorithms, we automatically insert at compile-time some compensation steps based on EFTs. More precisely, we introduce a tool named CoHD [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF] to parse C programs and generate a new C code with a compensated computation: floating-point operations ± and × are replaced by their respective error-free algorithms TWOSUM and TWOPRODUCT [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]Chap. 4]. Since using a generic solution (such as operator overloading) is not compatible with such transformations, CoHD generates a specific inline code.

We also present the SyHD software to perform transformations aiming to improve accuracy with an execution time requirement. SyHD synthesizes C source code for both accuracy and executiontime criteria. It uses CoHD transformations to improve accuracy, and transformation strategies to reduce how this accuracy improvement impacts the execution time. To improve execution time, we define two kinds of partial compensation strategies: some compensate the floatingpoint computations that generate the largest errors, and others do loop transformations, such as loop fission [START_REF] Appel | Modern Compiler Implementation: In ML[END_REF]. Then, SyHD synthesizes a new program by automatically analyzing the set of transformed programs for a given environment (defined by a target, some data, as well as some accuracy and execution-time constraints).

To demonstrate the efficiency of this approach, we compare the automatically transformed algorithms to existing compensated ones such as for floating-point summation [START_REF] Ogita | Accurate sum and dot product[END_REF] and polynomial evaluation [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF]. Our previous work CoHD [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF] led to automatically recovering the same level of accuracy and execution time. Compensation is the right choice to benefit from a high instruction level parallelism (ILP) compared to the solutions derived using fixed-length expansions such as double-double or quad-double [START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF][START_REF] Langlois | More instruction level parallelism explains the actual efficiency of compensated algorithms[END_REF]. The automatically transformed algorithms are shown to be very close, both in terms of accuracy and execution time, to the compensated algorithms we consider here. Now, our synthesizer SyHD applies transformation strategies to generate programs with accuracy and execution time trade-offs. We successfully generate programs with accuracy and execution-time trade-offs in a reasonable amount of time. We describe how this approach applies to a significant set of examples, ranging from recursive summation to polynomial evaluations and iterative refinement for linear system solving.

Related work

Several program transformation techniques or tools that improve the accuracy of floating-point computations have been proposed in the last decades. Nevertheless, most of them do not particularly take care of the execution time. One recent exception is Precimonius [START_REF] Rubio-González | Precimonius: Tuning assistant for floating-point precision[END_REF] that attempts to decrease the precision of intermediate operations to improve run-time and memory use. Extending the precision is a commonly used technique. Libraries that implement expansions, like the QD Library [START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF] or the Scilab toolbox [START_REF] Saito | Development of quadruple precision arithmetic toolbox qupat on scilab[END_REF], easily improve accuracy in simple cases. MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF] is an arbitrary-precision floating-point library with correct rounding, to accurately evaluate floatingpoint computations with the all necessary precision. These techniques introduce serious overheads for runtime and memory, as shown in [START_REF] Priest | On properties of floating point arithmetics: Numerical stability and the cost of accurate computations[END_REF] for example. Rewriting expressions. In [START_REF] Ioualalen | Synthesizing accurate floating-point formulas[END_REF], Ioualalen and Martel propose the Sardana tool. It performs a bounded exhaustive search for algebraically-equivalent programs for which a better accuracy bound could be proven statically. This approach is based on abstract interpretation to bound rounding errors using a reasonable over-approximation. A set of equivalent programs is generated over the real numbers, and the one with the smallest rounding error is chosen. Panchekha et al. [START_REF] Panchekha | Automatically improving accuracy for floating point expressions[END_REF] propose the Herbie tool, which improves the accuracy of expressions by randomly sampling inputs, localizing error, generating candidate rewrites, and merging rewrites with complementary effects. Herbie's results show that it can effectively discover transformations that substantially improve accuracy while imposing a median runtime overhead of 40% on their numerous selected case studies. A more aggressive solution, like the one used in GCC when using the -fast-math flag, speeds up floatingpoint programs by rewriting floating-point computations even if the numerical results may change. GCC gives no guarantee about the resulting accuracy but warns the user of potential numerical failures.

All these approaches tackle the question of accuracy improvement from different angles but share the tendency to combine this improvement with other criteria, especially with execution-time. Our tools can be used as a numerical analysis assistant but do not include the verification of numerical codes. Assisting and certifying numerical programs is for instance the goal of Gappa [START_REF] De Dinechin | Certifying floating-point implementations using gappa[END_REF], a tool based on the Coq proof assistant [START_REF]The Coq proof assistant reference manual[END_REF], or Fluctuat [START_REF] Goubault | Static analysis of numerical algorithms[END_REF], used to track floating-point errors with abstract interpretation.

Outline

This article is organized as follows. Section 2 provides the notations and some background material on floating-point arithmetic, EFTs, and accuracy improvement techniques like doubledouble arithmetic or compensation. The core of this article is Section 3 and Section 4, where we present our automatic code transformations to optimize the accuracy of floating-point computations with a small execution time overhead. Figure 1 gives an overview of this process. The four bottom modules compose the CoHD tool and the four top ones describe SyHD. Most of our contributions are implemented in the module "Transformation library," which is described in both Section 3 and Section 4. Section 4 also describes the modules "Strategies set research space," "Measurements & result exploration," and "Strategy selection." Other modules are classical compilation or code transformation blocks. In Section 5, some experimental results illustrate the practical interest of our approach compared to existing ones. Conclusions and perspectives are given in Section 6.

Input management

Strategies set research space

Measurements

& result exploration

PRELIMINARIES AND NOTATION

In this section we use classical notations to address IEEE floating-point arithmetic, basic methods to analyze the accuracy of floating-point computations, and EFTs of the basic operations ± and ×.

We also demonstrate how to exploit these EFTs within expansions and compensations.

IEEE floating-point arithmetic

In base 2 and precision p, IEEE floating-point numbers have the form

f = (-1) s • m • 2 e , where s ∈ {0, 1} is the sign, m = p-1 i=0 d i 2 -i = (d 0 .d 1 d 2 • • • d p-1
) 2 is the significand (with d i ∈ {0, 1}), and e is the exponent. The IEEE 754-2008 standard [START_REF][END_REF] defines such numbers for several formats, rounding modes, and the semantics of the five basic operations ±, ×, ÷, √ •.

Notation and assumptions. Throughout the paper, all computations are performed in the binary64 format, with the round-to-nearest mode. We assume that neither overflow nor underflow occurs during the computations. We use the following notations.

• F is the set of all normalized floating-point numbers. For example, the binary64 format, over 64 bits, includes p -1 = 52 bits for the fractional part of the significand, 11 bits for the exponent e, and 1 bit for the sign s. • f l(•) denotes the result of a floating-point computation where every operation inside the parenthesis is performed in precision p and with rounding to nearest. • Given x ∈ F, we write ulp(x) to denote the unit in the last place of x, defined by ulp(x) = 2 e • 2 1-p if x = 0, and by 0 otherwise. Let x = f l(x) for a real number x. We have |x -x| ulp( x)/2.

Accuracy analysis. One way to measure the accuracy of x = f l(x) is the number of significant bits # sig shared by x and x:

# sig ( x) = -log 2 |x -x|/|x| , x = 0. (1)

Error-free transformations

In this section, we briefly introduce the error-free transformations (EFT) previously discussed in [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF] (we refer to the latter for a more detailed description). EFTs provide lossless transformations of basic floating-point operations • ∈ {+, -, ×}. The practical interest of EFTs comes from FASTTWOSUM, TWOSUM (introduced by Dekker [1], Knuth [START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF]Chap. 4] and Møller [START_REF] Møller | Quasi double-precision in floating-point addition[END_REF]), and TWOPRODUCT [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] algorithms which exactly compute in floating-point arithmetic the error term of the sum and the product. Figure 2 defines diagrams for floating-point operations ± and ×, and for their EFTs. In what follows, we graphically represent transformation algorithms as basic computational blocks.

Double-double and compensated algorithms

We will now focus on two methods using these EFTs to double the accuracy: double-double expansions and compensations. We will then recall why compensated algorithms run faster than double-double algorithms (we refer to [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF] for a complete discussion on these two approaches). Double-double expansions. We consider here the Briggs, Kahan, and Bailey algorithms used in the QD library [START_REF] Hida | Algorithms for quad-double precision floating point arithmetic[END_REF]. Double-double arithmetic simulates computations with twice precision. Proofs are detailed in [START_REF] Lauter | Basic building blocks for a triple-double intermediate format[END_REF]. In practice, double-double algorithms can be simply derived by substituting the basic operations.

Compensated algorithms. As double-double algorithms, compensated algorithms can double the accuracy but they are not easy to derive. Compensated algorithms have been, up to now, defined on a case by case basis and by experts of rounding error analysis [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF][START_REF] Jiang | Accurate evaluation of a polynomial in Chebyshev form[END_REF][START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF][START_REF] Jiang | Accurate evaluation of a polynomial and its derivative in Bernstein form[END_REF][START_REF] Ogita | Accurate sum and dot product[END_REF]. For example, the compensated Algorithm SUM2 [START_REF] Ogita | Accurate sum and dot product[END_REF] returns a sum that is proven to be twice as accurate.

In this paper, we refer to a compensated algorithm by using the prefix "Comp" or the suffix "2," and to a double-double algorithm with the prefix or suffix "DD."

Double-double versus compensation. Double-double and compensated algorithms can be used to obtain similar accuracy. How do they compare in terms of computing time? A detailed analysis has been presented in [START_REF] Goossens | Computing time for summation algorithm: Less hazard and more scientific research[END_REF] concerning the classical summation algorithm. The authors show that SUM2 algorithm benefits from an ILP seven times higher than that of SUMDD. This fact is measurable in practice, and the compensated algorithms exploit this low-level parallelism much better than double-double ones. The example of HORNER's polynomial evaluation is detailed in [START_REF] Langlois | More instruction level parallelism explains the actual efficiency of compensated algorithms[END_REF]. The latter shows that the compensated HORNER's algorithm runs at least twice as fast as its double-double counterpart with the same output accuracy. This efficiency motivates us to automatically generate compensated algorithms.

COHD: AUTOMATIC ERROR COMPENSATION

We present solutions for improving accuracy thanks to code transformation, in order to benefit from the efficiency of compensation. This code transformation provides trade-off between accuracy and time compared to double-double expansions for example. More advanced trade-offs obtained with code synthesis will be discussed in Section 4.

CoHD is a source-to-source transformer written in OCaml and built as a compiler. The frontend reads input C files and comes from a previous development by Casse [START_REF] Casse | Frontc 3.4: an OCaml C parser and pretty-printer 2000[END_REF]. The middle-end implements some passes of optimization, from classical compiler passes such as operand renaming or three-address code conversion [START_REF] Appel | Modern Compiler Implementation: In ML[END_REF]Chap. 19]. It also implements one pass of floating-point error compensation we define in next Section 3.1. Then, the back-end translates the intermediate representation into C code.

Improving accuracy: methodology

Our code transformation automatically compensates programs following three steps.

1. Detection of sequences of floating-point computations. A sequence is the set S of dataflowdependent operations required to obtain one or several results. 2. For each sequence S i , computation and accumulation of the error terms. This is performed besides the original sequence by (a) replacing floating-point operations by the corresponding EFTS, and (b) accumulating error terms following Algorithms 1 and 2 given hereafter. At this stage, every floating-point number x ∈ S i becomes a compensated number, denoted x, δ x , where δ x ∈ F is the accumulated error term attached to the computed result x. 3. Closing of the sequences. Closing is the compensation step itself, so that close(S i ) means computing x ← f l(x + δ x ) for x being any result of S i .

Compensated operators for compensated numbers

Next Algorithms 1 and 2 automatically compensate the error of basic floating-point operations when inputs are compensated numbers. Two different sources of errors have to be considered.

First, the error generated by the elementary operation itself, which is computed with an EFT. This computation corresponds to δ ± and δ × in the first line of Algorithms 1 and 2. Second, the errors inherited from the two operands, denoted by δ a and δ b , are accumulated with the previous generated error. This accumulation corresponds to the second line of Algorithms 1 and 2. These inherited errors come from previous floating-point calculations. Operands with no inherited error are processed differently to minimize added compensations.

[s, δ ± ] = TWOSUM(a, b) Elementary operation replaced by the EFT

δ s ← f l((δ a + δ b ) + δ ± )
Accumulation of generated and inherited errors return s, δ s Algorithm 1: AC TWOSUM( a, δ a , b, δ b ), automatically compensated sum of two compensated numbers.

[s, δ × ] = TWOPRODUCT(a, b)
Elementary operation replaced by the EFT

δ s ← f l(((a × δ b ) + (b × δ a )) + δ × )
Accumulation of generated and inherited errors return s, δ s Algorithm 2: AC TWOPRODUCT( a, δ a , b, δ b ), automatically compensated product of two compensated numbers.

Application to computation sequence

In this section, we illustrate how our methodology transforms Listing 1 to Listing 2. First, CoHD detects the sequence a = b + c * d (Listing 1, line 4). This sequence is preliminarly converted in three-address form (Listing 2, lines 10 and 22). Then, computations are replaced by EFTs, and errors are accumulated following Algorithm 2 for the multiplication at line 10 (lines 11 to 21), and Algorithm 1 for the addition at line 22 (lines 23 to 27). Finally, the last step closes the sequence here at line 30 where the error term delta a is added to its corresponding variable. SyHD is a code synthesizer also written in OCaml. As described in Figure 1, it uses CoHD to perform multiple code transformations from some inputs being the C program file to optimize, optimization constraints, and data. SyHD explores a set of transformed programs generated with CoHD following several transformation strategies and finds the program with the best accuracy and execution time properties corresponding to the input parameters. Here, we describe how SyHD builds this transformed program set and how it selects the one with the best accuracy and time tradeoff.

c H = t -(t -c) ; c L = c -c H ; t = 134217729.0 * d ; d H = t -(t -d) ; d L = d -d H ; tmp L = c L * d L -((( tmp -c H * d H) -c L * d H) -c H * d L) ; / *

Automatic compensation

The automatic compensation process described in Section 3 improves accuracy with a minimal performance overhead compared to other techniques. Our SyHD tool here goes further to minimize execution-time overhead by partially applying the compensation of floating-point computations. Hence we reduce the performance overhead of the transformed program by sacrificing some Propagation of inherited errors return s, δ s Algorithm 3: AP TWOSUM( a, δ a , b, δ b ), sum with automatic error propagation of two compensated numbers.

s ← RN (a × b)
Elementary operation unchanged

δ s ← RN (a × δ b + b × δ a )
Propagation of inherited errors return s, δ s Algorithm 4: AP TWOPRODUCT( a, δ a , b, δ b ), product with automatic error propagation of two compensated numbers.

Again, operands with no inherited errors are processed to minimize added operations. Figure 4 shows such variants which can be obtained by removing the dashed or dotted lines. Here the elementary operations are not replaced by an EFT as for Algorithms 1 and 2 of Figure 3, so only inherited errors are processed (remove dashed or dotted lines for standard inputs).

One source of error remains now compared to previous Algorithms 1 and 2. Here, we only consider the inherited errors of the previous computations and ignore the generated error. Error propagation corresponds to the second lines of Algorithms 3 and 4 (the first lines are the unchanged elementary operations). This class of propagation is depicted in Figure 5b and it is denoted by p s which means "propagation with single closure." On the contrary, the closures performed at the end of each compensated part do not require the use of these additional algorithms. Computations not affected by compensation remain unchanged. We denote this case as p m , meaning "propagation with multiple closures." It is represented in Figure 5a.

Optimization strategies

In this section, we describe how programs are partially compensated. We propose two kinds of optimization strategies. A first pair of transformation operates on loops, while a third one operates on the floating-point computations which generate the largest errors. We use the minimal

i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 i = 7 i = 8 i = 9 i = 1 0
(a) Error propagation with multiple closures (denoted pm).

(b) Error propagation with single closure (denoted ps).

(c) Transformation SLT( , 1/3) with propagation pm.

(d) Transformation SLT( , 1/3) with propagation ps.

(e) Transformation ILT(⊥, 1, 4) with propagation pm.

(f) Transformation ILT( , 4, 5) with propagation ps.

(g) Transformation ILT(⊥, 2, 3) with propagation pm. We also define some basic procedures needed in Algorithms 5, 6, and 7 to manipulate program representation described in Listing 3. procedure has no effect if the propagation p m is in use, or applies automatic error propagation algorithms on floating-point sequences as explained in Section 4.1.

n = 1 n = 2 n = {3, 4} (h)
Simple loop transformation. SLT splits a loop into two loops to partially compensate one of them. Algorithm 5 defines the code transformation of this strategy. SLT depends on parameters (ρ, r), where r is the ratio of loop iterations to compensate, and ρ ∈ { , ⊥} is the position of these loop iterations. For example, SLT( , 1/2) compensates the first half iterations of the original loop. Figures 5c and5d Intermittent loop transformation. ILT splits a loop into 2k blocks in order to partially compensate one block every two blocks. The resulting program transformation consists in two loops nested into another one. Algorithm 6 describes this code transformation. ILT depends on three parameters (ρ, t, f ), where t is the number of loop iterations to compensate, f is the frequency of repetition (in number of loop iterations), and ρ ∈ { , ⊥} specifies the position of the loop iterations to compensate: the first t iterations of f if ρ = , the last t iterations otherwise. For example, ILT( , 1, 2) compensates the first iteration for each block of 2 iterations. Figures 5e,5f, and 5g illustrate this kind of transformations and an example is given in Listing 6.

Require: A loop statement id * min * max * body as defined in Section 4, and the strategy parameters: the number t of loop iterations to compensate, the frequency f (in iterations) of compensation, and the position ρ ∈ { , ⊥} Ensure: A loop statement replacing the required input statement m1 = expr stmt(j decrease order for all i ∈ {1, n} do compensate(get expr(pop(e))) get expr returns the expression statement of the computation which has generated the absolute error returned by pop end for return the input sequence with n compensated floating-point computations Algorithm 7: Accuracy oriented transformation (AOT).

+ t -1 > max ? j + t -1 : max) m2 = expr stmt(j + f -1 > max ? j + f -1 : max) if ρ == then replace t by f -t in m1, m2 l1 = loop stmt(id, j,

Strategies

Our loop transformation have to be sized with respect to the number of loop iterations. One additional parameter ν is introduced to be set by the user to adapt the transformations (ILT and AOT) to the program size. Figure 6 presents the default parameters for SLT, ILT and AOT. By default ν = 1. So, for example, if ILT operates on a loop performing one million iterations, then the default transformation splits the loop into blocks of size ranging from two iterations to five. By setting ν to 100 or 100,000, ILT splits the loops into blocks of sizes ranging from 200 to 500 or 200,000 to 500,000, and so it adapts the transformation to the number of loop iterations. Let us note that, where computational sequence begin after entering the loop, propagations p s and p m are equivalent for the SLT transformation when ρ = ⊥ (when appropriate the propagation is denoted by p m|s ). Then, the set E of strategies contains more than sixty different programs and

Code synthesis: accuracy and execution time trade-offs

We propose two ways to select the transformed program corresponding to accuracy and time constraints. R acc (E) selects the faster program among the most accurate ones and R exe (E) selects a less accurate program (compared to the one given by R acc ) to reduce the execution time. We introduce r sig. bits and r cycles respectively as the ratios of significant bits and cycles numbers and defined as: r * = automatically transformed program measure high-precision program measure .

(

) 2
We have: R acc (E) = min s∈{max s∈E r sig. bits (s)} r cycles (s), and,

R exe (E) = min s∈E α • r cycles (s) + β • (1 -r sig. bits (s)) α + β , ∀α, β > 0. (3)
R exe can be tuned to select a more or less accurate or fast program according to α and β. Moreover, to strictly satisfy the success of the selection, the transformed program must verifies the following points.

1. For each datum of a given dataset, the transformed program improves accuracy compared to the original program. 2. For the R acc criterion, the transformed program runtime must be smaller than the higherprecision program. 3. For the R exe criterion, the transformed program runtime must be smaller than the fully compensated program.

So, for a given program p and from: the set of strategies E, a set of data D, and an executiontime versus accuracy criterion R ∈ {R acc , R exe }, the program selection defined by Algorithm 8 returns a transformed program p ω satisfying the R criterion. Algorithm 8 introduces two procedures, get cycle number and get significant bit number, which return the number of cycles to execute a given program, and its accuracy expressed in the numbers of significant bits. Several solutions could implement these measures. We describe our preferred solutions in Section 5. Moreover, Algorithm 8 needs to compute some reference values such as n sig. bits the average number of significant bits obtained with a higher-precision version of the input program. In order to obtain the most efficient optimization strategy, data must be provided by the user. Although our tools do not support missing datasets, we could process our analysis on a substantial randomly chosen subset of the corresponding floating-point numbers set.

Algorithm 8 implies that the user has to preprocess the program to transform it by adding directives specifying which results have to be considered for accuracy improvement (typically, the result of a function). The user also has the option to specify the piece of code where time is measured.

EXPERIMENTAL RESULTS

We now describe the CoHD and SyHD tools that implement our code transformations. We apply them to several case studies described in Section 5.2 and chosen such that compensated versions are available to compare with and to validate the complete program transformation. We show in Section 5.3 that we automatically recover the accuracy and the performance of the existing compensated algorithms. We also add comparisons with the double-double versions. In Section 5.4 we explore the partial program transformation strategies in order to provide execution time versus accuracy tradeoffs. We show that we can generate transformed programs dealing with such Environment 2. AMD Athlon TM 64 X2 Dual Core Processor 4000+, Linux 3.2.0.51-generic-pae i686 athlon i386, gcc v4.6.3 with -O2 -mfpmath=sse -msse2, PAPI v5.1.1 and PERPI (pilp5 version).

In the following, Environment 1 is the default. Accuracy is measured as the number of significant bits # sig (see Relation ( 1)) in the floatingpoint significand. So, 53 is the maximum value we can expect from the binary64 format. For this study, we determine the accuracy with higher-precision reference values precomputed by MPFR within our tools.

A reliable measure of the execution time is more difficult to obtain. Such measurements are not always reproducible because of many unwanted side effects (operating system, executing programs. . . ). The most significant possible measures are provided here using two tools. The first one is the well known PAPI library (Performance Application Programming Interface) [START_REF] Mucci | PAPI: A portable interface to hardware performance counters[END_REF] which reads the hardware counters of cycles or instructions of an actual execution. The second software, PERPI [START_REF] Goossens | PerPI: A tool to measure instruction level parallelism[END_REF], measures the numbers of cycles and instructions of one ideal execution, that is, one execution by a machine with infinite resources. This measure is more related to a performance potential than to the actual measure provided by PAPI. Using both tools allow us to present confident and complementary results. Others measures are also presented in [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF].

Case studies

The first case study corresponds to existing compensated algorithms and is summarized with Table I and Table II. Table I gives data used for the summation of n values (case 1). Table II gives representative data used for the polynomial evaluation algorithms (cases 2, 3 and 4). Case 1. SUM2 for the recursive summation of n values [START_REF] Ogita | Accurate sum and dot product[END_REF]. Case 2. COMPHORNER [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF] and COMPHORNERDER [START_REF] Jiang | Accurate evaluation of the k-th derivative of a polynomial and its application[END_REF] for the Horner's evaluation of p H (x) = (x -0.75) 5 (x -1) 11 and its derivative. Case 3. COMPDECASTELJAU and COMPDECASTELJAUDER [START_REF] Jiang | Accurate evaluation of a polynomial and its derivative in Bernstein form[END_REF] for evaluating with de Casteljau's scheme p D (x) = (x -0.75) 7 (x -1) and its derivative, written in the Bernstein basis. Case 4. COMPCLENSHAWI and COMPCLENSHAWII [START_REF] Jiang | Accurate evaluation of a polynomial in Chebyshev form[END_REF] for evaluating with Clenshaw's scheme p C (x) = (x -0.75) 7 (x -1) 10 written in the Chebyshev basis.

Complete program transformation

We detail the complete program transformation of the Horner evaluation (case 2). Others cases will be summarized further.

Horner's polynomial evaluation case. We automatically compensate the Horner scheme and compare it with DDHORNER and COMPHORNER. The compensated algorithm and the data come from [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF]. Let p H (x) = (x -0.75) 5 (x -1) 11 be evaluated with HORNER's scheme, where x ∈ x 9 . Figure 7 shows the accuracy of this evaluation using HORNER (original), DDHORNER, COMPHORNER, and our automatically generated ACHORNER algorithm. In each case, we measure the number of significant bits # sig . The original HORNER algorithm's accuracy is low since the evaluation is performed in the neighborhood of multiple roots: most of the time, there is no significant bit. The other algorithms yield better accuracy. Our automatically generated algorithm exhibits the same accuracy behavior as DDHORNER and COMPHORNER, which are twice as accurate as the original algorithm. Table III shows the actual and ideal performances of the algorithms in terms of numbers of instructions, cycles, and instructions per cycle (IPC). The automatically generated algorithm has almost the same number of instructions and cycles as the compensated one. Moreover, Table III confirms that compensated algorithms benefit from more ILP than double-double ones. Even if the code of our generated algorithm is slightly different from the existing one [START_REF] Graillat | Algorithms for accurate, validated and fast polynomial evaluation[END_REF], it appears here to be as accurate and efficient. Moreover, the multi-criteria optimization introduced in Section 4 will produce quite different algorithms, by trading off accuracy against speed.

Others cases. We now summarize our results for the case studies of Table I and Table II. They have been chosen such that the original algorithm returns no significant digit at the working precision, while all of them are recovered by the twice more accurate ones.

Figure 8 presents the differences, in terms of number of significant bits, between the automatically compensated algorithms (AC) and the compensated ones (COMP) or the double-double ones (DD). For example, the eleventh case concerns HORNER's evaluation of p H for data x 1 . The difference between the numbers of significant bits of ACHORNER and COMPHORNER is zero. The AC algorithm is as accurate as the existing compensated one. The difference with the DD algorithm is of one bit. Most of the other results exhibit a similar good behavior. The slight differences of the last three data sets for the summation are due to the different effects of the sum length onto the compensated and double-double solutions: the accuracy bound of the former is quadratically affected by the length while being only linearly dependent in the double-double case. This appears only when the condition number and the sum length are large enough.

DECASTELJAUDER(p

D , x 2 ) DECASTELJAU(p D , x 1 ) DECASTELJAU(p D , x 2 ) CLENSHAWII(p C , x 1 ) CLENSHAWII(p C , x 2 ) CLENSHAWI(p C , x 1 ) CLENSHAWI(p C , x 2 )
HORNERDER(p

H , x 1 )
HORNERDER(p

H , x 2 )
HORNER(p

H , x 1 )
HORNER(p

H , x 2 ) SUM(d 1 ) SUM(d 2 ) SUM(d 3 ) SUM(d 4 ) SUM(d 5 ) SUM(d 6 )
AC-COMP AC-DD Figure 8. Differences of the numbers of significant bits # sig for the automatically compensated (AC) algorithm versus the existing compensated algorithms (COMP), and the double-double (DD) ones.

Figure 9 presents the performance ratios between AC algorithms and existing COMP or DD ones. Here, PAPI measures the actual mean of 10 6 executions. PERPI ideal measures are presented in [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF]. Figure 9a shows the cycle and instruction ratios. For example, the rightmost plot is SUM(d 8 ). We see that the instruction ratio AC/COMP = 1. Compensated SUM2 and our automatically generated one share the same number of instructions. The instruction ratio compared to DD is 0.8 which means that 20% fewer instructions have been generated while returning the same accuracy. The compensated algorithms, original and generated, present a similar number of cycles that remarkably represents 30% of the DD execution ones. Figure 9b shows the ratio of the number of instructions per cycle. We observe that AC algorithms have the same interesting properties as the original compensated ones. Measurements also confirm the interest of compensated algorithms, exhibiting a better ILP potential than DD algorithms. Finally, we note that the ILP potential is not fully exploited in our experimental environment [START_REF] Thévenoux | Automatic source-to-source error compensation of floating-point programs[END_REF].

Partial program transformations

We now consider the partial program transformations from Section 4.2 and 4.3. We begin with a detailed example of code synthesis using the HORNER algorithm. Then we summarize the results for the case studies of Table I andTable II. Horner's polynomial evaluation. Figure 10 summarizes the accuracy and performance improvements measured for all the transformed programs generated by the synthesis tool SyHD. 10 also displays how the accuracy and performance ratios (r sig. bits and r cycles ) vary according to the automatized transformations. We also add the original binary64 program and the completely compensated full-comp generated by CoHD. We see that the 62 strategies of partial transformation manage to provide an equivalent significant bits ratio (r sig. bits , closer to 1 the better) included from 0.6 to 0.7 and fail to achieve significant accuracy improvement compared to the original program (binary64), excepted for 5 SLTs strategies. The latter compensates for at least 70% (up to 90% for the best ratio) of the original program's computations, allowing to obtain from 77% to 88% of the accuracy of the completely compensated program full-comp. The cycles ratio (r cycles closer to 0 the better) distribution is quite different. We recover what we would intuitively expect: the more the computations are compensated, the more the number of cycles increases.

p m IL T( , 1, 2) , p s ILT ( , 1, 2), p m ILT(⊥, 1, 2), ps ILT (⊥, 1, 2), pm IL T( , 1, 3) , ps IL T ( , 1 , 3 ), pm IL T ( ⊥ , 1 , 3 ) , ps I L T ( ⊥ , 1 , 3 ) , pm I L T ( , 1 , 4 ) , p s I L T ( , 1 , 4 ) , p m I L T ( ⊥ , 1 , 4 ) , p s I L T ( ⊥ , 1 , 4 ) , p m I L T ( , 2 , 3 ) , p s I L T ( , 2 , 3 ) , p m I L T ( ⊥ , 2 , 3 ) , p s I L T ( ⊥ , 2 , 3 ) , p m IL T ( , 2 ,
Finally, the figure also points out the programs selected by the two trade-off selectors R acc and R exe (see relations ( 3) and ( 4)). Here SyHD default parameters apply (see Figure 6 and α = β for R exe ). R acc selects the fully compensated program (full-comp, depicted by •) because it is the only one which achieves the highest level of accuracy. This is due to the ill-condition of our data as explained for Table IV. R exe selects a program transformed following an SLT strategy: SLT(⊥, 9/10), p m|s (depicted by •). This transformation was selected because it allows the highest level of accuracy in the shortest execution time than the fully compensated one. As expected, an other program can be selected if R exe parameters was set differently.

Figure 11 gives the accuracy of the two selected programs from the Figure 10. It presents the number of significant bits when evaluating p(x) = (x -0.75) 5 (x -1) 11 with the Horner's scheme and x in the x 3 dataset. Figure 11a compares the accuracy of the original and transformed programs when the R acc selector is desired. Full accuracy is recovered for approximately 289 cycles which is more than 3 times faster than the double-double version. Figure 11b shows the accuracy behavior when the R exe selector is required. By sacrificing accuracy (from 0 to 12 significant bits as shown by the "output" line in Figure 11b), we can save at least 13% of cycles for this polynomial evaluation (which now takes 251 cycles).

This example illustrates the potential of our approach. From Figure 11, we can consider the generation of programs presenting accuracy between the same level of accuracy as the original and the full-comp (completely compensated) programs. These programs could be generated by tweaking the selector parameters by requesting some different accuracy and execution time trade-offs. 10) and 11b (SLT(⊥, 9/10) with p m|s , also depicted by • in Figure 10). For all these cases, the R acc selector returns the completely compensated algorithm as the transformed solution. This case is explained by the ill-condition of the tested case: no transformation except the fully compensated one yields a fully accurate compensation. Less ill-conditioned data would be transformed differently by R acc . The R exe selector only yields partially compensated programs with SLT transformation mainly because of the ill-condition of the data. Nevertheless, other partial transformations could give an appropriate result if R exe was set differently or if data was less ill-conditioned.

Other cases.

Table IV also displays some failures, denoted by † or ‡ when selecting a program with R exe . The following notes explain them. † This failure is raised because the accuracy requirement is not reached. As defined in the requirements of the code synthesis success (Section 4), the transformed program must improve the accuracy for each desired datum. This requirement is very strong and we could modify it (to make it customizable) to accept that the optimization as no effect on accuracy for some data. We could also customize the R exe selector in order to accept a less accurate trade-off by setting α = β/2 for example. By doing so, the case denoted by a returns the transformed program obtained with the strategy SLT(⊥, 7/10), p m|s . ‡ This failure is raised because of a lack of performance. DECASTELJAU's algorithms are tiny codes (a loop of 8 iterations containing a few floating-point operations), and this explains that every partially compensated program doesn't perform better than the fully compensated one.

In such cases, failure can be avoided by returning this latter.

Finally, we claim that this code synthesis is processed in a reasonable amount of time. It takes from about 1 minute (Environment 1) to 1.5 minute (Environment 2) for the polynomial evaluation cases. In these cases, the code synthesis has to explore more than 60 different programs, over 256 data to measure accuracy, and between 100,000 and 1,000,000 executions to measure performance. For the summation cases, with larger dataset and longer execution time, the code synthesis can take up to 2 minutes (Environment 2) to perform the program selection.

Iterative refinement for linear system solving

Finally, we consider the classical iterative refinement technique used for improving the accuracy of the solution of a linear system Ax = b. Algorithm 9 performs such computation [START_REF] Demmel | Error bounds from extra-precise iterative refinement[END_REF]. When computing the residual in twice the working precision (line 4), iterative refinement returns the fully accurate solution of not too ill conditioned linear systems [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]Chap. 12]. The BLAS [START_REF] Blackford | An updated set of Basic Linear Algebra Subprograms (BLAS)[END_REF] (Basic Linear Algebra Subroutines) routines are convenient here since they are ready to use and welltuned for performance. The XBLAS library [START_REF] Li | Design, implementation and testing of extended and mixed precision blas[END_REF] (extended and mixed precision BLAS) relies on BLAS and also provides some routines performing double precision inner computation, e.g. iterative refinement. The goal of this experiment is to compare our transformed code computing the residual to the one of the XBLAS.

Require: An n × n matrix A and an n × 1 vector b Ensure: A solution vector x (i) approximating x in Ax = b 1: Solve Ax Solve Ax (i+1) = r (i)

6:

Update x (i+1) = x (i) -x (i+1)

7:

i = i + 1 8: until x (i) is "accurate enough" 9: return x (i) Algorithm 9: Classical iterative refinement.

Figure 12 shows the iteration number and the accuracy of Algorithm 9 when the matrix A is orthog [START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF] or gfpp(50) for three variants of the algorithm: the residual is computed with BLAS, XBLAS, or the AC generated code. The elements of the vector b are uniformly distributed in [0, 1]. The stopping condition (line 8) is satisfied when ||x (i+1) -x (i) ||/||x (i+1) || 4u, where u = 2 -53 . The data of this experiment are obtained with the Matrix Computation Toolbox [START_REF] Higham | Matrix Computation Toolbox[END_REF] and chosen similarly to [35, pp. 239-240]. Figure 12 exhibits that the XBLAS implementation and AC generated code share both the same performance and accuracy. Table V presents the average cycle measure of 20 executions of Algorithm 9. It confirms that the XBLAS and AC codes present equivalent performance and perform refinement with a similar number of cycles. In order to beat XBLAS performance, we will focus on a better code transformation from CoHD as discussed in Section 6. In this article we discussed the automated transformation of programs using floating-point arithmetic. We propose a new method to compensate automatically the floating-point errors of the computations, and to improve the accuracy without greatly impacting execution time. The automatic transformation produces some compensated algorithms which are as accurate and efficient as the ones derived by hand on a case-by-case basis. Moreover, we propose a partial compensation to reduce further the execution time overhead by trading off performance and accuracy. Trade-offs are ensured by code synthesis and multi-criteria program optimization. We developed strategies to generate partially transformed programs, as well as a method to find the best transformation satisfying an execution time or accuracy constraints. The efficiency of our approach has been illustrated in various case studies. We now need to validate this approach and the tools by testing them on real and more sophisticated programs. To achieve this, we have to support floating-point divisions, square-roots, and the elementary functions. Possible solutions for improving our tools are numerous. First, we would like to add more transformation strategies in order to generate new partially compensated program patterns, or other optimization constraints such as code size or energy consumption. Then, we aim to develop a dynamic method for Algorithm 8 to find a program that satisfies the optimization constraints, for example, by using a genetic algorithm. We would also like to yield accuracy by adding an unbounded compensation step applying the step proposed in this paper several times, or by using expansions and compensation. Finally, in order to facilitate these perspectives, we plan to implement the CoHD transformations into the GCC compiler.

Figure 1 .

 1 Figure 1. CoHD (bottom) and SyHD (top) tools detailed in Section 3 and Section 4.

Figure 2 .

 2 Figure 2. Diagrams for basic floating-point operations (a), (b) and EFT algorithms (c), (d), and (e).

Figure 3

 3 Figure3shows such variants which can be obtained by removing the dashed or dotted lines.

Figure 3 .

 3 Figure 3. Diagrams for Algorithms 1 and 2 for the automatic compensation of the sum (a), and the product (b). The dashed or dotted lines are removed when a or b is a standard floating-point number.

Listing 1 :

 1 Original code computing the sequence a = b + c × d. double foo() { double a, b, c, d ; [...] a = b + c * d ; / * computation sequence * / 5 [...] return a ;

Listing 2 :

 2 Transformed code computing the sequence a = b + c × d with error compensation.

  step 2b: accumulation of TwoProduct error term result and inherited errors from c and d * / delta tmp = tmp L ; / * second part of the sequence detected at step 1 * / a = b + tmp ; / * step 2a: adding 5 flops with TwoSum(b,tmp) * / t = ab ; a L = (b -(at)) + (tmpt) / * step 2b: accumulation of TwoSum error term result and inherited errors from b and tmp * / delta a = a L + delta tmp; [...] / * step 3: close sequence * / return a + delta a ; } 4. SYHD: CODE SYNTHESIS TO OPTIMIZE ACCURACY AND TIME Section 3 explains how to generate a fully compensated program. Here we introduce strategies of partial transformation that yield programs with different patterns of transformation. Hence we aim to select between them the best trade-off for accuracy and runtime constraints defined by the user.

Figure 4 .

 4 Figure 4. Diagrams for Algorithms 3 and 4 for the automatic error propagation of the sum (a), and the product (b).Here the elementary operations are not replaced by an EFT as for Algorithms 1 and 2 of Figure3, so only inherited errors are processed (remove dashed or dotted lines for standard inputs).

  A possible transformation AOT(4) with propagation ps.

Figure 5 .

 5 Figure 5. Representation of transformed program traces (sequential execution). Figures 5a and 5b define the basics of the representation: compensated parts are shaded and error propagation is represented with a sub bold line. Figures 5c to 5h show several examples of program transformation.

10 L

 10 . TH ÉVENOUX ET AL.

Algorithm 5 :

 5 Simple loop transformation (SLT).

SLT r = δ 1 10Figure 6 .

 16 Figure 6. Default parameters of strategy set E for code synthesis. Parameter ν can be chosen by the user to adapt the transformations to the program size (ν = 1 by default).

Figure 7 .

 7 Figure 7. Number of significant bits # sig (y axis) when evaluating p H (x) = (x -0.75) 5 (x -1)11 , wherex ∈ x 9 (x axis) for HORNER, DDHORNER, COMPHORNER, and ACHORNER.

Figure 9 .

 9 Figure 9. Performance ratios between automatically compensated algorithms (AC) and existing compensated (COMP) or double-double (DD) ones (mean of 10 6 values measured with PAPI).

  4/10), ps SLT ( , 4/1 0), pm SL T( ⊥ , 4/ 10 ), pm|s S L T ( , 5 / 1 0 ), ps S L T ( , 5 / 1 0 ) , pmS L T ( ⊥ , 5 / 1 0 ) , pm |s S L T ( , 6 / 1 0 ) , p s S L T ( , 6 / 1 0 ) , p m S L T ( ⊥ , 6 / 1 0 ) , p m | s S L T ( , 7 / 1 0 ) , p s S L T ( , 7 / 1 0 ) , p m S L T ( ⊥ , 7 / 1 0 ) , p m | s S L T ( , 8 / 1 0 ) , p s S L T ( , 8 / 1 0 ) , p m S L T ( ⊥ ,8 / 1 0 ) , p m |s S L T ( , 9 / 1 0 ), p s SL T( , 9/ 10 ), p m SLT (⊥, 9/1 0), p m|s • ILT( , 1, 5), ps ILT (

Figure 11 .

 11 Figure 11. Accuracy (significant bits numbers) for (a) Racc and (b) Rexe selection of the transformed HORNER's evaluation for the x 3 data (Table II, x axis).

(1) = b 2 : i = 1 3: repeat 4 :

 24 Compute residual r(i) = Ax (i) -b 5:

Figure 12 .

 12 Figure 12. 20 executions (x axis) of Algorithm 9 where A is matrix orthog(25) or gfpp(50). The y axis shows the iteration number until termination (left) and the number of lost digits (right).

  When partially compensating a program, new alternatives occur. Indeed, after a compensated part, the accumulated errors can be propagated through the next floating-point computations instead of performing a closure. For addition and multiplication, this propagation is achieved by Algorithms 3 and 4 given below. Elementary operation unchanged δ s ← RN (δ a + δ b )

Copyright c 0000 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (0000) Prepared using cpeauth.cls DOI: 10.1002/cpe accuracy. s ← RN (a + b)

  illustrate it and an example of such transformation is given in Listing 5.Require: A loop statement id * min * max * body as defined in Section 4, and the strategy parameters: the ratio r of iterations to compensate, and the position ρ ∈ { , ⊥} Ensure: A sequence statement replacing the required input statement if ρ == then e = f loor((maxmin) × r) l1 = loop stmt(id, min, e + min, compensate(body)) l2 = loop stmt(id, e + min + 1, max, no compensate(body)) else

	end if
	return sequence(l1, l2)

e = f loor((max -min) × (1 -r)) l1 = loop stmt(id, min, e + min, no compensate(body)) l2 = loop stmt(id, e + min + 1, max, compensate(body))

  Listing 4. Recursive summation in C. int main(int argc, char * * argv) { int i, n=atoi(argv[1]); double s, tab[n]; FILE * file; file = fopen (argv[2], "rb"); Accuracy oriented transformation. AOT selects the floating-point operations that generate the largest errors. Errors can be measured through different techniques. We propose to consider the absolute errors of the transformed code results compared to some reference values computed in higher precision. The principle of this transformation is explained in Algorithm 7 and Figure 5h illustrates it. Require: A program p to transform, a set D of dataset, and the strategy parameter: the number n of expression statements to compensate Ensure: The transformed input program which compensates the n expressions of p generating the maximum absolute errors compared to its higher-precision counterpart for all d ∈ D do execute the program p and its counterpart in higher precision with the data d P d ← all the intermediate results of the floating-point computations of p P d ← all the intermediate results of the floating-point computations of p in higher precision (reference results) for all ( i, i) ∈ ( P d , P d ) do

		fread(&tab, sizeof(double), n, file);
		fclose (file);
		s = tab[0];
		for(i=1; i<n; i++)
		s = s + tab[i];
	}	printf("%a\n", s); e d i ← E abs ( i, i) end for return 0; end for vector e contain all the absolute errors of all the intermediate floating-point the absolute error i -i computations of p for each dataset in D
		e ← mean(e)	component-wise
	Listing 5. Listing 4 after transformation e ← sort(e)	Listing 6. Listing 4, after transformation ILT( , 1, 2)
		SLT( , 1/2) with propagation pm.	with propagation ps.
	int main(int argc, char * * argv) {	int main(int argc, char * * argv) {
		int i, n = atoi(argv[1]); double s, tab[n] ;	int i, n = atoi(argv[1]); double s, tab[n] ;
		int ac end;	int ac i, ac end2, ac end3;
		double s dt, ac tmp1, ac tmp1 dt,	double s dt, ac tmp1, ac tmp1 dt,
		s twosum, s ts1, s ts2;	s twosum, s ts1, s ts2;
		FILE s = tab[0]; s dt = 0.0; ac end = (int) n * 0.5; for(i = 1; i < ac end; i++) { ac tmp1 = s; ac tmp1 dt = s dt;	FILE for(i = ac i; i < ac end2; i++) {
		s = ac tmp1 + tab[i];	ac tmp1 = s;
		s ts1 = s -tab[i];	ac tmp1 dt = s dt;
		s ts2 = s -s ts1;	s = ac tmp1 + tab[i];
		s twosum = ( ac tmp1 -s ts1)	s ts1 = s -tab[i];
		+ (tab[i] -s ts2);	s ts2 = s -s ts1;	20
		s dt = s twosum + ac tmp1 dt;	s twosum = ( ac tmp1 -s ts1)
		}	+ (tab[i] -s ts2);
		s += s dt;	s dt = s twosum + ac tmp1 dt;
		s dt = 0.0;	}
		for(i = ac end; i < n; i++) {	for(i = ac end2; i < ac end3; i++) {
		s = s + tab[i];	s = s + tab[i];
		}	s dt = s dt;
			}
		printf("%a\n", s + s dt);	}
	}	return 0; l2 = loop stmt(id, m1 + 1, m2, no compensate(body)) m1, compensate(body)) printf("%a\n", s + s dt); return 0;	30
		else
		l1 = loop stmt(id, j, m1, no compensate(body))
		l2 = loop stmt(id, m1 + 1, m2, compensate(body))
		end if
		e = expr stmt(j = j + f )
		return loop stmt(j, min, max, sequence(l1, l2, e))

* file; file = fopen(argv

[START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF]

, "rb"); fread(&tab, sizeof( double ), n, file); fclose(file); * file; file = fopen(argv

[START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF]

, "rb"); fread(&tab, sizeof( double ), n, file); fclose(file); s = tab[0]; s dt = 0.0; for(ac i = 1; ac i < n; ac i = ac i + 2) { ac end2 = ac i + 1 < n ? ac i + 1 : n; ac end3 = ac i + 2 < n ? ac i + 2 : n; }

  large set of possible partially compensated programs. Section 4.4 describes how to select the best program into E to match a required accuracy and time trade-off.

Copyright c 0000 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (0000) Prepared using cpeauth.cls DOI: 10.1002/cpe AUTOMATIC ERROR COMPENSATION OF FLOATING-POINT PROGRAMS 13 hence covers a

Table I .

 I Case studies and data for the summation (case 1).

		Table II. Case studies and data for polynomial evalua-
		tion with HORNER, CLENSHAW, and DECASTELJAU
				(cases 2, 3, 4).
	Summation Data # values condition number d 1 32 × 10 4 10 8 d 2 32 × 10 5 10 8 d 3 32 × 10 6 10 8 d 4 32 × 10 4 10 16 d 5 32 × 10 5 10 16 d 6 32 × 10 6 10 16 d 7 10 5 irrelevant d 8 10 6 irrelevant	Polynomial evaluations Data # x range x 1 256 {0.85 : 0.95} (uniform dist.) x 2 256 {1.05 : 1.15} (uniform dist.) x 3 256 {0.35 : 0.45} (uniform dist.) x 4 256 {0.6 : 0.7} (uniform dist.) x 5 256 {1.8 : 1.9} (uniform dist.) x 6 256 {1.2 : 1.3} (uniform dist.) x 7 256 {0.73 : 0.74} (uniform dist.) x 8 256 {0.755 : 0.8} (uniform dist.) x 9 512 {0.68 : 1.15} (uniform dist.)
		x	1	irrelevant

Table III .

 III Performance measurements of the algorithms COMPHORNER, DDHORNER, and ACHORNER. Real values (PAPI) are the means of 10 6 measures. Ideal values are obtained with PERPI.

			PAPI			PERPI
		instructions cycles IPC instructions cycles IPC
	COMPHORNER	532	277	1.99	566	62	9.12
	DDHORNER	658	920	0.72	676	325	2.08
	ACHORNER	553	303	1.82	581	77	7.54
	AC/COMP	1.04	1.09	0.95	1.01	1.24	0.82
	AC/DD	0.84	0.33	2.55	0.85	0.23	3.62

  Results exploration in the case of HORNER's algorithm with data x 3 (r cycles closer to 0 the better, r sig. bits closer to 1 the better, • denotes the program returned with Racc, and • the one returned with Rexe). Figure 10 corresponds to the HORNER's algorithm running on data x 3 (256 values in {0.35 : 0.45}). Our transformation strategies (SLT, ILT, AOT) generate here 62 different programs.Figure

	4 ) , p s	IL T ( , 2 , 4 ), p m 2, IL T( ⊥ , 4) , p s	ILT (⊥, 2, 4), p m	ILT( , 3, 4), ps	ILT ( , 3, 4), pm	IL T( ⊥ , 3, 4) , ps IL T (⊥ , pm , ), 4 3 A 1 ) O ( T	• p m l f u l -o c b i n y 4 6 r a	r cycles r sig. bits
	Figure 10. Copyright c 0000 John Wiley & Sons, Ltd.					Concurrency Computat.: Pract. Exper. (0000)
	Prepared using cpeauth.cls							DOI: 10.1002/cpe

  Table IV summarizes the code synthesis results obtained for all the case studies of Table I and Table II. Default SyHD parameters apply except for SUM(d 1 ) and SUM(d 2 ) where ν is respectively set to 1,000 and 10,000. The bold results correspond to the ones depicted in Figures 11a (full-comp, also depicted by • in Figure

Table IV .

 IV Results summarize of code synthesis on all the case studies (TableIand TableII). The bold results correspond to the ones depicted in Figures 11a (full-comp, also depicted by • in Figure10) and 11b (SLT(⊥, 9/10) with p m|s , also depicted by • in Figure10).

	Algorithms and data	Optimization criteria	Optimization criteria
		Environment 1	Environment 2
		Racc	Rexe	Racc	Rexe
	SUM(d 1 )	full-comp	fail †	full-comp	fail †
	SUM(d 2 )	full-comp	fail †	full-comp	fail †
	HORNER(p H , x 3 )	full-comp SLT(⊥, 9/10), p m|s full-comp	SLT( , 9/10), pm
	HORNER(p H , x 4 )	full-comp	SLT( , 9/10), ps	full-comp	SLT( , 9/10), pm
	HORNERDER(p H , x 3 )	full-comp	fail †	full-comp	fail †
	HORNERDER(p H , x 4 )	full-comp	SLT( , 9/10), ps	full-comp	SLT( , 9/10), ps
	CLENSHAWI(p C , x 3 )	full-comp	SLT(⊥, 8/10), p m|s	full-comp SLT(⊥, 7/10), p m|s
	CLENSHAWI(p C , x 5 )	full-comp	SLT( , 8/10), ps	full-comp SLT(⊥, 8/10), p m|s
	CLENSHAWII(p C , x 3 )	full-comp	SLT(⊥, 8/10), p m|s	full-comp SLT(⊥, 8/10), p m|s
	CLENSHAWII(p C , x 5 )	full-comp	SLT( , 8/10), ps	full-comp SLT(⊥, 8/10), p m|s
	DECASTELJAU(p D , x 7 )	full-comp	fail ‡	full-comp	fail ‡
	DECASTELJAU(p D , x 8 )	full-comp	fail ‡	full-comp	fail ‡
	DECASTELJAUDER(p D , x 7 ) full-comp	fail ‡	full-comp	fail ‡
	DECASTELJAUDER(p D , x 8 ) full-comp	fail ‡	full-comp	fail ‡

Table V .

 V Minimal, average and maximal kilocycles (PAPI) of 20 executions of Algorithm 9. XBLAS AC min. ave. max. min. ave. max.

	orthog(25) 483 552	790	496 568	816
	gfpp(50)	109 114	180	112 116	186
	6. CONCLUSIONS AND PERSPECTIVES

Copyright c 0000 John Wiley & Sons, Ltd. Concurrency Computat.: Pract. Exper. (0000) Prepared using cpeauth.cls DOI: 10.1002/cpe

constraints, and this in a reasonable amount of time. Test cases include summation, polynomial and derivative evaluations with HORNER, CLENSHAW, and DECASTELJAU algorithms for which compensated versions have been published. A last test case presented in Section 5.5 is the iterative refinement algorithm applied to several ill-conditioned linear systems.

We recall that compensated and double-double algorithms are respectively prefixed with "Comp" and "DD." Similarly, we use the "AC" prefix for those generated by our automatic compensation method. We start to introduce how we perform time and accuracy measurements as required by Algorithm 8.

Experimental environments and measure methods.

We perform accuracy and execution time measurements to compare programs automatically generated from our method to hand-coded compensation programs and double-double algorithms. All measurements are done within the following experimental environments.


Sébastien Lefebvre
email: sebastien.lefebvre@vetagro-sup.fr

Benoit Rannou

Stéphane Besse

Marlène Damin-Pernik

Etienne Benoit

Virginie Lattard

Preliminary study on the efficiency of vitamin K on hypocoagulable rats. Introduction

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction:

Vitamin K are essential molecules involved in the post-translational activation of some proteins called protein vitamin K dependent (PVKD). The best know PVKD are the four involved in the coagulation cascade: factor II, VII, IX and X. In mammals, vitamin K has two origins: foods and a recycling system led by an enzyme called Vitamin K epoxide reductase complex subunit 1 (VKORC1). The vitamin K requirement is dependent of the efficiency of this enzyme [1]. The aim of our work is to block the recycling mechanism in order to evaluate the efficiency of vitamin K 1 to balance the recycle system and so to better understand the vitamin K 1 action.

Results:

The delay between difethialone and vitamin K1 administration have been chosen to maximize the diminution of the vitamin K dependent clotting factor activities without involving death [2]. The results of the regeneration of vitamin K dependent clotting factor activities are presented on the figure 2. Regeneration evolutions fitted with a sigmoidal model (Factor II R²=0.88; VII R²=0.84; IX R²=0.78; X R²=0.93). After 360 minutes, factor II and factor IX were significantly lower (P<0.05) than their basal activity level before difethialone administration (set at 100%). For all factors, a lag time was observed during the first hour after the vitamin K administration. Then, activities increased during about one hour before to stagnate. The prothrombin time was within the normal range from 120 minutes (figure 3). The lag time was expected but not the different levels of factor activities. The difference between factors might be explained by different affinities in the post-transcriptional modification. The stop of increasing might be explained by lack of vitamin K. This arrest of regeneration was not described in the case of regeneration by recovery of recycling mechanism [3]. So, unrecycled vitamin K might have a different activity or affinity which leads to a new equilibrium.

Conclusion:

Our experiment pinpoints that the regeneration of clotting

  factor activities is more complex to understand than expected. The study of the clotting factors' regeneration might disclose some new regulation mechanisms. By means of data provided by this study we are able to design new experiments, in order to describe properly the regeneration mechanisms. References: 1. Markussen MD, Heiberg A-C, Nielsen R, Leirs H. Vitamin K requirement in Danish anticoagulant-resistant Norway rats (Rattus norvegicus). Pest Manag Sci. 2003;59(8):913-20. 2. Lefebvre S, Rannou B, Besse S, Benoit E, Lattard V. Origin of the gender differences of the natural resistance to antivitamin K anticoagulants in rats. Toxicology. 17 févr 2016;344-346:34-41. 3. Kerins GM, MacNicoll AD. Comparison of the half-lives and regeneration rates of blood clotting factors II, VII, and X in anticoagulant-resistant and susceptible Norway rats (Rattus norvegicus Berk.). Comp Biochem Physiol C Pharmacol Toxicol Endocrinol. mars 1999;122(3):307-16.

Figure 3 :Figure 1 :Figure 2 :

 312 Figure 3: Evolution prothrombin time after administration of Difethialone then, 48h after, administration of Vitamin K