Clothilde Jeangoudoux
email: clothilde.jeangoudoux@safrangroup.com

Christoph Lauter
email: christoph.lauter@lip6.fr

A Correctly Rounded Mixed-Radix Fused-Multiply-Add

The IEEE 754-2008 Standard governs Floating-Point Arithmetic in all types of Computer Systems. The Standard provides for two radices, 2 and 10. It specifies conversion operations between these radices, but does not allow floatingpoint formats of different radices to be mixed in computational operations. In contrast, the Standard does provide for mixing formats of one radix in one operation.

In order to enhance the Standard and make it closed under all basic computational operations, we propose an algorithm for a correctly rounded mixed-radix Fused-Multiply-and-Add (FMA). Our algorithm takes any combination of IEEE754 binary64 and decimal64 numbers in argument and provides a result in IEEE754 binary64 and decimal64, rounded according to any for the five IEEE754 rounding modes.

Our implementation does not require any dynamic memory allocation; its runtime can be bounded statically. We compare our implementation to a basic mixed-radix FMA implementation based on the GMP Multiple Precision library.

I. INTRODUCTION

When developing software using floating-point arithmetic, the user is presented with a choice between two radices: binary floating-point arithmetic (radix-2) or decimal floatingpoint arithmetic (radix-10). The IEEE 754-2008 Standard for Floating-Point Arithmetic [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] defines the format and operations for both radices, and states that a programming environment may conform to this standard in one radix or in both. Quite a few studies have been undertaken on decimal floatingpoint arithmetic in an effort to provide accurate and efficient decimal operations using binary formats [START_REF] Harrison | Decimal Transcendentals via Binary[END_REF] and conversion operations [START_REF] Cornea | A Software Implementation of the IEEE 754R Decimal Floating-Point Arithmetic Using the Binary Encoding Format[END_REF]. Those algorithms provide faithfully (or comparably accurate) rounded results. We would like to go one step further and perform correctly rounded "mixed-radix" floatingpoint arithmetic, that is, arithmetic operations with operands and result represented in a binary or a decimal floating-point format.

To perform mixed-radix operations, the user may be tempted to combine binary and decimal floating-point arithmetic such as in the following C code: When trying to compile this code e.g. with gcc 6.3.0, the compiler refuses to translate the code, emitting an error message which states that floating-point types of different radices cannot be mixed. The user, who might have no other choice than using and mixing two scientific libraries, one of which produces results in binary floating-point and the other produces decimal ones, might hence be inclined to force the compiler to translate this piece of code, e.g. by adding explicit casts (conversions). This attempt may also be misleading, as these conversions induce errors, which may be insignificant but may also accumulate, amplify and become important.

Extending our example code as follows showcases this dangerous [START_REF] Wong | Recent Catastrophic Accidents: Investigating How Software was Responsible[END_REF] behavior:

1 double a=2000.0, c=-2.0, d; 2 _Decimal64 b=0.001D; This code compiles, but does not set the correct result to d, i.e. 2000×0.001-2 = 0, but 4.16333634 •10 -17 . This error is due to the conversion of the decimal64 number b into a binary floating-point number.

All these reasons point toward the need for a study and conception of a specific mixed-radix framework. However, as of now, the applicability of mixed-radix operations has sparsely been studied and very few of them are available to the user. Previous work investigated the feasibility of exact comparison between floating-point formats of different radices [START_REF] Brisebarre | Comparison Between Binary and Decimal Floating-Point Numbers[END_REF], as a way to enrich the current floating-point environment and make numerical software safer. We would like to extend this approach to the basic computational operations defined by the IEEE754 Standard: +, -, ×, /, √ and FMA. The FMA operation, added in 2008, computes a × b + c with one, single, correct rounding. We wish to provide algorithms for a mixedradix implementation of the computational operations, i.e. algorithms that take floating-point numbers of both radicestwo and ten-supported by the IEEE754-2008 Standard and compute correctly rounded floating-point results, also in one of these two radices.

A naive approach would be to perform the FMA computation in one of those two radices and perform a rounding of the result into the output format. Suppose we have

a = 2 -55 • 3242591731706757, b = 2 -52 • 4548635623644201, c = 2 -52 • 5902958103587057.
All these three numbers are IEEE754 binary64 numbers. The binary64 result of the FMA operation is then

• (a × b + c) = 2 -52
• 5902958103587057 that rounds to the decimal64 10 -17 • 9090999999999999, where in a correctly rounded mixed-radix arithmetic the decimal64 result would have been 10 -17 • 9091000000000000. It is easy to construct an example for which the mixed-radix result of an operation will be hard to round. Yet we have no simple systematic way to seek inputs for which this will be the case, hence our need for a study and specific implementation of a correctly rounded mixed-radix FMA. Our work is based on a previous study [START_REF] Kupriianova | Towards a Modern Floating-Point Environment[END_REF], where an attempt to compute such hard-to-round cases was made, without actually succeeding.

There are several binary and decimal floating-point formats available, as defined in the IEEE754-2008 such as binary32, binary64 and decimal64 or decimal128, and in principle even mixed-radix heterogeneous operations like a binary64 plus decimal128 gives binary32-addition, might eventually to be considered. To get started, it is however easier to look at one binary and one decimal format, with comparable accuracy, i.e. precisions that are comparable provided appropriate conversion. Hence to perform those algorithms, we choose the IEEE754 binary64 format, which provides k = 53 bits of binary precision, which is equivalent to a basic relative error of u = 2 -53 , and the IEEE754 decimal64 that provides k = 16 digits of decimal precision, which is equivalent to a basic relative error of u = 1/2 × 10 -15 = 2 -50.82... . We will suppose that the IEEE754 decimal formats, such as the decimal64 format we use, are encoded in the IEEE754 binary encoding (BID) [START_REF] Cornea | A Software Implementation of the IEEE 754R Decimal Floating-Point Arithmetic Using the Binary Encoding Format[END_REF], i.e. that access to a binary representation of the decimal significand is straightforward. Given this setting, let us admit the following:

• Binary, resp. decimal, operands are denoted 2 E • m, resp.

10 F • n, with a bounded signed binary exponent E ∈ Z, resp. decimal exponent F ∈ Z, and m ∈ Z, a signed integer significand satisfying

2 k-1 ≤ |m| ≤ 2 k -1,
where k is the binary precision, resp. n ∈ Z with 1 ≤ |n| ≤ 10 l -1 and l the decimal precision. • The result of the computation of the mixed-radix FMA is either a binary output 2 G • p, with G ∈ Z appropriately bounded and p ∈ Z a signed integer significand bounded, or a decimal output 10 H • q with H ∈ Z appropriately bounded and q ∈ Z a signed integer significand bounded. The IEEE 754 Standard defines the notion of "quantum" for decimal operations [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. However, whenever an operation's input is binary (floating-point or integer), the quantum is defined in an ad hoc manner. We shall hence not consider it.

Throughout this paper, we will first highlight the challenges due to the computation of correctly rounded mixedradix operations (Section II). We will then describe the internal mechanisms involved in the computation of a correctly rounded mixed-radix FMA (Section III), and then discuss the implementation and test strategy developed to measure the efficiency of our approach (Section IV).

II. CHALLENGES AND SETTINGS

Provided a correctly rounded mixed-radix FMA such as • (a × b + c), we can easily compute the correctly rounded mixed-radix multiplication by setting the c operand to zero, and the addition/subtraction by setting either a or b to one.

Furthermore both correctly rounded mixed-radix division and square root can be determined with a single mixed-radix FMA operation, provided it allows for slightly larger precision: given a binary precision k, resp. a decimal precision l, we call midpoint, noted f , a number that is exactly halfway between two consecutive floating-point numbers.

The rounding of the result of the mixed-radix division • k x y of binary precision k, resp. decimal precision l, can be determined with the help of the closest midpoint of the real result x y [START_REF] De Dinechin | On Ziv's Rounding Test[END_REF]. This midpoint f can be represented with a floating-point number of slightly larger precision k , resp. l . Thus the rounding decision boils down to decide whether

• x y < f ⇔ 0 < y × f -x, • x y > f ⇔ 0 > y × f -x, • or x y = f ⇔ 0 = y × f -x.
This implies that identifying the rounding of the division x y and thus computing the correctly rounded mixed-radix division operation in the output format of precision k, resp. l, is equivalent to compute a mixed-radix FMA operation y × f -x that allows for operands of slightly larger precision k , resp. l , and compare its result to zero.

Similarly, the computation of the correctly rounded mixedradix square root of precision k, resp. l, boils down to determine where the result falls according to a midpoint f . Thus to determine if √ x < f we can compute 0 < f × f -x with a mixed-radix FMA operation that allows for operands of slightly larger precision k , resp. l .

For the binary64 format, we can represent all binary floating-point numbers, normal and subnormal, and their midpoints by increasing the binary precision from k = 53 to k = 55. The modification of the precision of the binary64 format has an impact on the bounds of its significand and exponent, thus it can be defined as follows

2 E • m; with 2 54 ≤ |m| < 2 55 ; -1130 ≤ E ≤ 969; m, E ∈ Z. ( 1
)
Since 10 is divisible by 2, the representation of the midpoints is also possible in the decimal64 format, by adding only one binary bit of precision, for the midpoint between two decimal floating-point numbers can be represented in the same binade instead of decade. This means that we can represent the significand n with the same 55-bit precision as the binary64 format earlier and adapt the exponent, by factorizing 10 F as 2 J • 5 K , with K = F . We can then adapt the decimal64 input representation such as

10 F • n = 2 J • 5 K • r; with 2 54 ≤ |r| < 2 55 ; r ∈ Z -452 ≤ J ≤ 385; -421 ≤ K ≤ 385; J, K ∈ Z. (2)
There is however not one mixed-radix FMA operation, but several operations, depending on the combination of formats given as inputs and output. As a matter of fact, the FMA operation handles as inputs the three numbers a, b and c, and returns one output, all four numbers being representable in two different formats: binary64 or decimal64. We can naively consider that we will need to implement 2 4 = 16 versions of the mixed-radix FMA. Yet as we are studying mixed-radix arithmetic, we can easily discard the all binary and all decimal FMA operations that will not benefit from our approach. Furthermore we will consider the decimal64-to-binary64multiplication and the binary64-to-decimal64-multiplication to be equivalent, whether we want a binary or decimal output. Hence this excludes four more versions of the mixed-radix FMA operation as redundancies. In the end, we will only consider ten different versions of the mixed-radix FMA operation.

The real challenge of mixed-radix operations lies into the computation of the addition or subtraction inside the FMA operation. As floating-point subtraction is just floating-point addition with a change in sign of the subtrahend, we are going to continue the discussion only for "addition", subsuming subtraction. Typically, for a floating-point addition, two cases can be distinguished: the case when the two input operands are sufficiently close in terms of exponent that the subtraction of their aligned significands will lead to cancellation, and the other case when they are far enough in terms of magnitude so that the result's order of magnitude is essentially the one of the bigger of the operands. The first of these two cases is commonly called near path, the latter far path [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF].

Uniform-radix correctly rounded floating-point addition is based, for the near path, on the observation that when cancellation occurs, the floating-point operation becomes naturally exact, i.e. does not require any rounding [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. This case corresponds to the case when Sterbenz' lemma can be applied [START_REF] Sterbenz | Floating-Point Computation, ser. Prentice-Hall series in automatic computation[END_REF]. For mixed-radix addition, there is no reason why this cancellation case should become exact: e.g. when both input operands are in decimal floating-point arithmetic, both 10 and -9.9 are representable. There sum of course is 0.1 which clearly is not representable in binary floatingpoint arithmetic. A decimal plus decimal to binary mixedradix addition can hence not suppose that no rounding occurs in near path cases. Furthermore, when the two operands to a mixed-radix addition are not in the radix, some sort of conversion of the binary operand to the decimal radix or vice versa will be needed. While cancellation in the near-path makes the operation become exact for uniform-radix cases, the same cancellation in the mixed-radix case will amplify this conversion error:

a + b • (1 + ε) = (a + b) • 1 + b a + b • ε , |a + b| |b| .
And even if we find ways to compute lower bounds on this catastrophic cancellation results of mixed-radix addition on its near path, i.e. if we find ways to make the conversion sufficiently accurate that some accuracy is left once the cancellation amplifies the error, the result of this (approximate) mixedradix addition may itself be close to a rounding boundary of the output radix, preventing us to return the correctly rounded result with certainty. On the far path, mixed-radix addition is no less difficult, as absorption can occur and lead to the same rounding issues. Even though we can construct exemplary cases where correct rounding of mixed-radix addition is hard to perform, we have no simple, systematic way to look for inputs to mixedradix addition for which rounding is hard in round-to-nearest or a directed rounding mode. We therefore propose another approach to implementing a correctly rounded mixed-radix FMA, without precomputing such lower bounds.

III. MIXED-RADIX FMA ALGORITHM

The input of the mixed-radix FMA algorithm are either binary64 or decimal64 with the possibility of representing the midpoints by slightly increasing the precision. Those two formats are defined by the equations ( 1) and ( 2), and can be unified with a single representation of the mixed-radix unified format, with pessimist bounds on the exponents such as, given a an input of the FMA operation, we have:

a = 2 Na • 5 Pa • t a ; with 2 54 ≤ |t a | < 2 55 ; t a ∈ Z -1130 ≤ N a ≤ 969; -421 ≤ P a ≤ 385; N a , P a ∈ Z. (3)
The algorithm that performs a fast and correctly rounded FMA defined in Figure 1 can be split in several parts. The first step is the decomposition of the binary64 and decimal64 inputs into the explicit format of sign, exponents and significand. Throughout this algorithm, we will analyze the bounds of those variables since we want to ensure that they can be stored on machine words: for the exponent on a signed 32-bit integer machine word and for the significand on one or several unsigned 64-bit integer machine words. With this internal format, we compute an error-free mixed-radix multiplication. Then we perform a test to decide which addition/subtraction algorithm to use, described in more detail later on. According to the result of this test, we convert the operands into an internal binary format and perform either a far-path addition/subtraction, or a near-path subtraction. We will then test if the final rounding in the output radix is possible with this fast-computed result. If the result is too close to a midpoint f , we cannot decide what is the correct floating-point value that should take the result according to the current rounding mode [START_REF] De Dinechin | On Ziv's Rounding Test[END_REF]. In that case we will have an extra computation step, called the recovery phase, to get those last bits of precision needed to decide how to round the result.

In the following Section, we will describe further the mechanisms of the mixed-radix FMA algorithm. We will first describe the far-path and near-path addition/subtraction algorithms of the fast computation, and then explain how work the rounding test and the recovery phase.

A. Multiplication and Ratio Test

The first operation performed by the FMA is pretty simple: the multiplication can be computed exactly in the mixedradix unified format. According the definition of the mixedradix unified format of the input operands a and b given in equation ( 3), we can compute the result of the multiplication such as

ψ = a × b = 2 L • 5 M • s; with 2 109 ≤ |s| < 2 110 ; s ∈ Z -2261 ≤ L ≤ 1938; -842 ≤ M ≤ 770; L, M ∈ Z.
(4) We then want to perform the addition/subtraction between ψ and c, the third input operand of the FMA that is represented in the format given in equation [START_REF] Cornea | A Software Implementation of the IEEE 754R Decimal Floating-Point Arithmetic Using the Binary Encoding Format[END_REF]. Previously, in Section II, we outlined two cases for the floating-point addition, according to the proximity of the two operands in terms of exponents. Those two cases are the far-path and near-path addition algorithms.

Deciding whether the operands are close or far enough to perform either one of those algorithm can determined by Require: a, b and c, binary64 or decimal64 numbers

Ensure: R = 2 G • p in binary or R = 10 H • q in decimal 1: ψ = a × b 2: if it is an "addition" or ψ c / ∈ [ 1 2 ; 2] then 3: T 1 ← 2 γ1 • v 1 = ψ • (1 + ε T 1 F P ) 4: T 2 ← 2 γ2 • v 2 = c • (1 + ε T 2 F P ) 5: φ = (T 1 ± T 2 )(1 + ε φ F P ) ("far-path") 6: else 7: T 1 ← 2 Γ1 • w 1 = ψ • (1 + ε T 1 N P ) 8: T 2 ← 2 Γ2 • w 2 = c • (1 + ε T 2 N P ) 9: φ = (T 1 -T 2 )(1 + ε φ N P ) ("near-path") 10: end if 11: ρ ← 2 C • g = φ or ρ ← 10 H • 2 -10 • q = φ • (1 + ε ρ D ) 12: if ρ = (a × b + c)(1 + ε ρ ) can round correctly then 13:
return R ← ρ correctly rounded into output format 14: else 15:

Compute integer rounding boundary significand f such as

2 C • 2 10 • 1/2 • f or 10 H • 2 -10 • 2 10 • 1/2 • f 16: Let α ← 2 L-Zmin • 5 M -Fmin • (a × b) + 2 Nc-Zmin • 5 Pc-Fmin • c -2 F T -Zmin • 5 F F -Fmin • f ) with α ∈ Z 17:
Correct ρ using f and the sign of α 18:

return R ← ρ correctly rounded into output format 19: end if ; 2], the operands are far enough to perform the faster far-path algorithm, otherwise the more accurate near-path algorithm is executed. This test ensures that when the result states that the far-path algorithm should be used the ratio is surely outside of

[ 1 2 ; 2].
Given the definitions of ψ and c given in equations ( 3) and (4), we want to ensure that either we have

2 L-Nc •5 M -Pc • s tc < 1 2 or 2 L-Nc • 5 M -Pc • s tc > 2.
Thus certifying that the ratio ψ c is clearly outside of [ 1 2 ; 2] is equivalent to evaluating the following inequation

(L -N c ) + (M -P c ) • log 2 (5) + 1 + 56 < -1 or (L -N c ) + (M -P c ) • log 2 (5) + 54 > 1 (5)
Given the precomputed value log 2 (5), A•log 2 (5) with A ∈ Z in a small domain can be computed correctly as A • log 2 (5) • 2 48 •2 -48 which does not overflow on a 64-bit signed integer variable [START_REF] Brisebarre | Comparison Between Binary and Decimal Floating-Point Numbers[END_REF]. Hence we can easily compute on a 64-bit signed integer the two values given in equation ( 6) and perform ratio test that gives us two cases.

B. When no cancellation occurs (far-path)

1) Far-path algorithm: When the ratio test states that the result is clearly outside of [ 1 2 ; 2], no cancellation can occur during the the addition/subtraction, hence we can use the far-path addition algorithm. For simplicity reasons, this algorithms will take binary floating-point inputs and return a binary floatingpoint result. As the operands ψ and c are represented in the mixed-radix unified format with precisions of respectively 110bit and 55-bit precisions. Thus we first want to perform a conversion of those operands into a 64-bit binary format. This means that we want to represent the significand v stored on a 64-bit unsigned integer such as 2 63 ≤ |v| ≤ 2 64 -1, and adapt the binary exponent γ, as described by the lines 3 and 4 of algorithm 1. The next step being the computation of the correctly rounded in the current rounding mode binary addition/subtraction and leave the rounding of the result into the output format for further considerations.

As the operands ψ and c do not have the same significand precision, we might first think that we need two conversion algorithms from unified mixed-radix format to a 64-bit binary format. Yet it can easily be shown that we can first unify those two operands into a 64-bit mixed-radix format such as

2 λ • 5 µ • ω; with 2 63 ≤ |ω| ≤ 2 64 -1; ω ∈ Z -2215 ≤ λ ≤ 1284; -842 ≤ µ ≤ 770; λ, µ ∈ Z, (6)
an then apply a unique conversion algorithm from this mixedradix format to a binary variable T represented with a 64-bit unsigned integer significand such as

T = 2 γ • v; with 2 63 ≤ |v| ≤ 2 64 -1; v ∈ Z -4171 ≤ γ ≤ 3772; γ ∈ Z. (7)
Applying this conversion algorithm on ψ and c gives us the two variables T 1 = 2 γ1 • v 1 and T 2 = 2 γ2 • v 2 defined as in equation ( 7) as inputs for the far-path addition/subtraction algorithm, described in Figure 2.

The overall philosophy of this algorithm is that when neither of the operands are zero, we can order T 1 and T 2 according to their exponents such as γ 1 ≥ γ 2 . Absorption occurs when γ 1 ≥ γ 2 + 64, in that case the result of the far-path addition is

2 γ 1 • v 1 .
Otherwise, to perform the addition of the two 64-bit unsigned integer exponents v 1 and v 2 , we pose the intermediate 128-bit unsigned integer variables V 1 and V 2 and align the exponents by shifting v 1 to the left such asV 1 = 2 γ 1 -γ 2 • v 1 . After performing the addition or subtraction on a 128-bit unsigned integer variable V = V 1 ± V 2 , the significand is normalized with the function lzc(V ) that counts the number of leading zeros and truncated to a 64-bit unsigned integer variable g. The exponent C is set accordingly.

The result of the far path addition algorithm is then returned in the following binary format

2 C • g; with 2 63 ≤ |g| ≤ 2 64 -1; g ∈ Z -4174 ≤ C ≤ 3836; C ∈ Z. (8)
2) Far-path error analysis: We will analyze the global relative error made during the computation of the far-path addition/subtraction algorithm, from the conversion of the inputs ψ and c until the acquisition of the result φ = 2 C • g. We shall denote this global error ε φ F P .

The first operation producing an error occurs during the computation of the equation ( 6). In the case of ψ, as the precision decreased from 110 bits to 64 bits, an error arises,

2 λ ψ • 5 µ ψ • ω ψ = 2 L • 5 M • s • (1 + ε ψ )
, and it can be bounded by |ε ψ | ≤ 2 -63 . Then, given an input in the format described in equation ( 6), we compute the exponent γ and significand Require: The global relative error attached to the computation of the 64-bit precision binary variables T 1 and T 2 is

v as γ = λ + µ • log 2 (5) ; v = 2 -64 • ω • τ 0 (µ) , with the precomputed table τ 0 (µ) = 5 µ •2 -µ•log 2 (5) +63 ,
T 1 = 2 γ1 • v 1 and T 2 = 2 γ2 • v 2 Ensure: φ = (T 1 ± T 2 )(1 + ε φ F P ) = 2 C • g 1: Order T 1 and T 2 so that γ 1 ≥ γ 2 2: if γ 1 ≥ γ 2 + 64 then 3: g ← v 1 4: C ← γ 1 5: else 6: V 1 ← 2 γ 1 -γ 2 • v 1 7: V 2 ← v 2 8: V ← V 1 ± V 2 9: σ ← lzc(V ) 10: π ← 2 σ • V 11: g ← 2 -64 • π 12: C ← γ 2 + 64 -σ 13: end if 14: return φ = 2 C • g
2 γ •v = 2 λ •5 µ •ω• (1+ε v )•(1+ε τ0(µ)
). The relative error of v is a combination of the truncation error and the error attached to the computation of the table τ 0 (µ) such as ε v is bounded by

|ε v | ≤ 2 -62 .
The global relative error for this part of the mixed-radix FMA algorithm includes the conversion error ε ψ , the error on the precomputed table ε τ0(µ) and the rounding error ε v . With those three errors combined, we deduce a pessimist bound on the global relative error |ε T F P | ≤ 2 -60. 5 . Given T 1 and T 2 the inputs of the far-path addition, we want to bound the relative error ε φ F P as φ = (T 1 ± T 2 )(1 + ε φ F P ).

The algorithm can compute the result in two different ways, firstly if

γ 1 ≥ γ 2 + 64, the result is φ = 2 γ 1 • v 1 • (1 + ε φ F P ) with |ε φ F P | ≤ 2 -63
. Secondly, if we are in the case when γ 1 < γ 2 +64, hence γ 2 ≤ γ 1 ≤ γ 2 +63 because γ 1 , γ 2 ∈ Z, we can deduce the bounds on V 1 and V 2 defined in the algorithm 2 such as 2 63 ≤ V 1 < 2 127 and 2 63 ≤ V 2 < 2 64 . We want to compute bounds on V = V 1 ± V 2 . Either it is an addition and we have 2 64 ≤ V 1 + V 2 < 2 127 + 2 64 < 2 128 , or it is a subtraction, and we can show, with a reductio ad absurdum, that 49 . Following the steps of the algorithm, we then bound σ according to its definition such as 0 ≤ σ ≤ 66; thus π is between 2 127 ≤ π ≤ 2 128 and finally we have

2 61 ≤ 1 2 •(1-2 -59.49 )•2 63 ≤ V 1 +V 2 < 2 127 +2 64 < 2 128 , with the fact that 1 -2 -59.49 ≤ 1+ε T 1 F P 1+ε T 2 F P ≤ 1 + 2 -59.
g = 2 -64 •π = 2 -64 • π • ε φ F P . In the end we have |ε φ F P | ≤ 2 -63 .

C. When cancellation does occur (near-path)

When the ratio test states that the result could be inside [ 1 2 ; 2], cancellation may occur during the computation of the subtraction. Hence the use of the near-path subtraction.

A solution for this near-path algorithm is to use the same method as the far-path one and perform the subtraction on binary operands. To do so, we need to determine the number of bits of precision necessary for the binary inputs to compute the near-path subtraction in the worst case of cancellation and avoid loosing all accuracy.

1) Worst case of cancellation: The worst case occurs when, without being zero, the subtraction between ψ and c such as Proof. By application of the algorithm found in [START_REF] Brisebarre | Comparison Between Binary and Decimal Floating-Point Numbers[END_REF].

2 L • 5 M • s -2 Nc • 5 Pc • t c is relatively small. This is equivalent to 2 L • 5 M • s ≈ 2 Nc • 5 Pc • t, thus s t ≈ 2 Nc-L • 5 Pc-M ,
2) Near-path algorithm: Subsequently, we deduce that to compute the near-path subtraction with at least 60 bits of precision, we need to represent the binary significands with at least 178 + 60 = 228 bits of precision. This means that we can represent the significands of the operands of the nearpath subtraction on a 256-bit unsigned integer, i.e. four 64-bit unsigned integers machine words.

Hence we first compute the conversion from the mixed-radix unified format variables φ and c to 256-bit precision variables, with which we perform the near-path subtraction.

As in the far-path algorithm, we want to perform only one conversion into the binary internal format, therefore adapt the significand and exponent of c by increasing to the 110-bit precision of ψ as in equation ( 4). Then perform the conversion from the 110-bit precision mixed-radix format to the 256-bit precision binary format such as

T = 2 Γ • w; with 2 255 ≤ |w| ≤ 2 256 -1; w ∈ Z -4364 ≤ Γ ≤ 3578; Γ ∈ Z. (9)
The near-path subtraction algorithm is quite similar to the far-path except that the alignment of the significands yields to an error. We first order the two operands T 1 = 2 Γ1 • w 1 and T 2 = 2 Γ2 • w 2 such as Γ 1 > Γ 2 . With the result T 1 and T 2 , we align the significands so that z 2 = 2 Γ 2 -Γ 1 • w 2 ; z 1 = w 1 . After ordering z 1 and z 2 such as z 1 > z 2 , we can compute the subtraction d = z 1 -z 2 , perform the count of leading zeros in d with lzc(d) such as we can normalize the significand and exponent of the result g = d • 2 l • 2 -192 ; and C = Γ 1 -l + 192. The result of the near-path subtraction takes the same form as the result of the far-path addition defined in equation ( 8), except for the bounds of the exponent that are -2296 ≤ C ≤ 3773.

3) Near-path error analysis: The first error that occurs during the computation of the near-path algorithm comes from the conversion from the mixed-radix unified format to the 256bit precision binary format. We can denote this error ε T N P and define it such as

2 Γ • w = 2 L • 5 M • s • (1 + ε T N P ), computed as Γ = L -18 + M • log 2 (5) + σ with 0 ≤ σ ≤ 1; σ ∈ N and w = 2 -129 • τ 0-3 (M ) • s • 2 18 • 2 σ .
The table τ 0-3 (M ) is the 256-bit precision version of the table τ 0 (µ) used previously. We can show that τ 0 (M ) = 2 -192 • 2 255-M •log 2 (5) • 5 M = 2 -192 • τ 0-3 (M ) . We can then bound the relative error induced by the computation of the table such as |ε τ0-3 | ≤ 2 -255 , and the global relative error of the conversion algorithm such as |ε

T N P )| ≤ 2 -251 .
The result φ of the near-path binary subtraction is subject to an error such as φ = (T 1 -T 2 )(1+ε φ N P ). A first error appears during the computation of If we want a binary output for the FMA operation, we can directly perform the rounding test that informs us whether the result can be correctly rounded into a binary64 floating-point number. If we want a decimal64 output, we must convert this 64-bit binary number into a decimal one such as

z 2 = (2 Γ 2 -Γ 1 • w 2 ) • (1 + ε z2 ). By noticing that -2 ≤ Γ 2 -Γ 1 ≤ 0 we can bound ε z2 such as |ε z2 | ≤ 2 -253 . This implies that g = (d•2 l •2 -192 )•(1+ε φ N P ) with 2 255 ≤ d • 2 l ≤ 2 256 -1.
10 H • 2 -10 • q with -1253 ≤ H ≤ 1159; H, q ∈ Z 10 15 • 2 10 ≤ q ≤ (10 16 -1) • 2 10 < 2 64 -1. ( 10
)
This algorithm computes the exponent

H = (C + 63) • log 10 (2) -15 + µ(C, g) where µ(C, g) ∈ N; µ(C, g) = C • log 10 (2) + log 10 (g) -(C + 63) • log 10 (2)
. We can show that 0 ≤ µ(C, g) ≤ 1, implying that µ(C, g) ∈ {0, 1}, by observing that µ(C, g) is an increasing function in g. Thus for a given C, there exists g * (C) ∈ Z;

2 63 ≤ g * (C) < 2 64 , µ(C, g) = 0 if g ≤ g * (C) 1 if g > g * (C). (11)
We can precompute such a g * (C) with a table of 64-bit unsigned integer words.

Once the exponent H is set, we can compute the significand q with a bounded precomputed table t(H) = 10 -H • 2 H•log 2 (10) +127 and multiply it by g. The value is then normalized and truncated properly, conforming to the definition, in equation 10, of this operation's output.

We skipped the detailed error analysis of this algorithm for the sake of conciseness but the full pencil-and-paper proof will soon be available under the form of a research report. The relative error of this conversion is bounded such as 10

H • 2 -10 • q = (2 C • g) • (1 + ε ρ D ) with |ε ρ D | ≤ 2 -59.82 .
2) Rounding test: The purpose of this test is to determine whether the 64-bit precision output of the addition/subtraction algorithm can be rounded into a a binary64, resp. decimal64, with any sign and rounding mode.

In the binary64 case, we have as input

2 C • g = (ψ ± c)(1 + ε * B ) with ε * B the
C • g such as 2 F T • 5 F F • f = 2 C • 2 10 • 1 2 • f and 2 54 ≤ f ≤ 2 55 -1,
and step into the recovery phase.

In the decimal64 case, we have as input

2 H •2 -10 •q = (ψ ± c)(1+ε * D )
-2 -9 • q 2 9 | ≥ β D .
If this condition is not satisfied, rounding to a 16-digit precision decimal64 number is not possible. Therefore we compute f , the binary rounding boundary closest to 2 H • q such as 2

F T • 5 F F • f = 10 H • 2 -10 • 2 10 • 1 2 • f and 2 • 10 15 ≤ f ≤ 2 • (10 16 -1)
, and proceed to the recovery phase.

3) Recovery phase: If the rounding test fails, we try to round a result that is too close to the rounding border f to be able to decide towards which nearest floating-point number it must be rounded. Similarly to the computation of this floating-point rounding border's value in the previous step, the recovery phase boils down to computing the sign of

α = (-1) sa • 2 L-Zmin • 5 M -Fmin • s + (-1) sb • 2 Nc-Zmin • 5 Pc-Fmin • t c -2 F T -Zmin • 5 F F -Fmin • f with bounds on the exponents such as 0 ≤ F T -Z min ≤ 2162; 0 ≤ F F -F min ≤ 786; 0 ≤ L -Z min ≤ 2282; 0 ≤ M -F min ≤ 785; 0 ≤ N c -Z min ≤ 2342; 0 ≤ P c -F min ≤ 786.
Given those bounds, we know that we need to compute positive or zero integer powers of 5, i.e. 5 k , with 0 ≤ k ≤ 786. Furthermore, we know that every 5 k holds on log 2 (5 786 ) = 1826 bits, hence 29 words of 64 bits. Finally, we determine that each term of α holds on respectively on 4218, 4223 and 4043 bits, so whatever their sign is, they hold on a maximum of 4225 bits. So to compute exactly α and decide its sign, we need 67 words of 64 bits, leaving 63 "free" bits.

We have three cases: if α = 0 then the result is exactly the rounding border f , and will be rounded in the output format correctly by our final rounding algorithm according to the current rounding mode. If α < 0, then the result is below the rounding border f and will round to the lower nearest floatingpoint number in the output format. By analogy, if α > 0 then the result shall round to the higher nearest floating-point number in the output format. Computing this recovery phase is therefore very costly, but fortunately this part of the FMA algorithm will only be computed sparsely.

4) Final rounding to the output format: Final rounding to the supported output formats, binary64 and decimal64, takes either (-1) s • 2 C • g, 2 63 ≤ g < 2 64 , or (-1) s • 10 H • 2 -10 • q, 10 15 • 2 10 ≤ g < 10 16 • 2 10 , as input and rounds them to a binary64 or decimal64 floating-point number respectively, according to the current IEEE754 rounding mode. Binary rounding is done in the IEEE754 rounding mode for binary floating-point arithmetic, decimal rounding in the decimal mode. This final rounding also sets the inexact, underflow or overflow flags, when the rounding is respectively inexact, underflowed (tiny and inexact) or overflowed.

The difficulty in implementing this final operation is twofold: first, the rounding must address all special cases, such as subnormal rounding, normal rounding, overflow, decimal denormalization at the least exponent etc. Second, the rounding must be performed with the appropriate rounding mode, access to which is difficult. It also must set the flags. We overcome both challenges by exactly computing binary64, resp. deci-mal64, numbers that, provided to a IEEE754 floating-point operation, round to the same value as the input data.

For binary64 rounding, we decompose (-1) s 11 • g . For decimal64 rounding, we represent (-1) s • 10 H • 2 -10 • q exactly as an IEEE754 decimal128 floating-point value and perform an IEEE754 decimal128 to decimal64 format conversion, which rounds correctly and sets the flags. This way of performing final decimal rounding is known to be sub-optimal in terms of performance but has the advantage of being easy to implement. We consider changing this part of our algorithm and perform decimal64 rounding manually, while detecting the appropriate rounding mode by executing the avatar decimal64 instruction t = (-1) s • 10 15 + 1 resp. t = (-1) s • 10 15 + 2 added to δ, with δ ∈ {-3/4, -1/2, -1/4, 0, 1/4, 1/2, 3/4}, choosing t, t and δ depending on the parity of a truncated significand and the value of the rest that got truncated off.

• 2 C • g. Typically C = C + 11, g = g • 2 -11 , A + B = C, r = g -2

IV. IMPLEMENTATION AND TEST RESULT

A. Reference implementations

In order to perform software testing on our mixed-radix FMA implementation, we designed two other implementations of a mixed-radix FMA: one, based on the GNU Multiple Precision Library (GMP) [START_REF] Granlund | the GMP development team, GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF], with the same C calling interface as our optimized FMA implementation, and one other, written in Sollya [START_REF] Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF], with an ad hoc text-based interface. While designing these two reference implementations, we set the following design goals: the GMP-based implementation should be as close as possible to the result of a mixed-radix FMA designed in a limited timeframe, reusing existing software library. It should be reasonably fast, but easy to design.

We addressed this design goal by converting all input operands, be they IEEE754 binary64 or decimal64 numbers, into GMP rational numbers, i.e. fractions of (long) integers. We then performed the multiplication and addition steps of the FMA in GMP rational numbers. GMP can perform these operations on rational numbers without any error, computing exact long products, sums and gcds as required, while allocating memory dynamically. Our reference implementation then rounds the rational number, call it p q , by first determining the output exponent for a binary64, resp. decimal64, number, call it E, and then "rounding" the rational number r s = 2 -E • p q (resp.

r s = 10 -E • p q ) to an integer significand by long division of the long integer r by s, adapting the division rest's sign as required by the rounding mode. As the input/output's exponents may grow relatively large, GMP will manipulate integers of size of about a couple of thousand bits in the worst case.

In the next Section IV-B, we give an overview of the testing of our optimized FMA implementation we performed, comparing its results to the GMP-based reference implementation, and also compare timing results. While performing software development and intermediate testing of both our optimized implementation and the GMP-based one, we found certain input operands where the two algorithms did not agree. In order to have a way to decide where to look out for issues, we designed yet another implementation of a mixed-radix FMA, in the interpreted numerical language Sollya [START_REF] Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. In this tool, numerical expressions can be written exactly, and evaluated to any precision without being hurt by spurious roundings. It is rather easy to give definitions of the output exponent and significand of a mixed-radix FMA in terms of the input exponents and integer significands. Particularities due to subnormal rounding, underflow and overflow can be addressed through minimum/maximum computations.

We used our Sollya reference implementations to generate "gold" result test vectors, with which we then verified both our optimized FMA implementation and the GMP-based reference. After our software testing campaign, all our three implementations agreed. While we spent quite some time on testing the implementation of our mixed-radix FMA algorithm against reference implementations, we would like to insist that we consider our algorithm to be correct only because we were able to prove its correctness, the detailed pencil-and-paper proof made available under the form of a research report.

B. Test and Timing Results

We tested our mixed-radix FMA implementation in all its 14 different operand-result-format variants on a system based on an Intel i7-7500U quad-core processor, clocked at maximally 2.7GHz, running Debian/GNU Linux 4.9.0-5 in x86-64 mode. We compiled our code using gcc version 6.3.0 with optimization level 3 and setting the -march=native flag. On our system, this results in AVX instructions being used, including hardware binary64 FMA. The IEEE754-2008 decimal types and instructions come from our gcc's libgcc. The GMP-based reference implementation was based on GMP version 6.1.2. Sollya was at git commit 12bf2006c.

We tested each of the 14 variants on an appropriate test vector file containing at least 115000 test points, covering all possible cases of signs, signed zeros, infinities, NaNs, subnormal and normal numbers, as well as all possible binary and decimal rounding modes 1 . We found all our 14 mixedradix FMA implementations to be correct with respect to their numerical outputs for all the outlined cases.

We also performed testing of the implementations concerning the correct setting of the IEEE754-2008 flags, as they are inexact, overflow, underflow, divide-by-zero and invalid [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. We found all flags to be set correctly, i.e. our implementations correctly raised the appropriate IEEE754-2008 when the operation actually was inexact (i.e. implied rounding), overflowed, underflowed or was invalid (as for multiplication of zero by infinity). The implementation also was checked for correctly leaving the flag unchanged if it did not encounter the suited IEEE754-defined condition to raise the flag. One issue was found about this last check for the mixed-radix FMA implementations for certain decimal exact results, the inexact flag was raised spuriously. We traced this undue raising of the inexact flag back to the implementation of the IEEE754 decimal128 to decimal64 format conversion that is contained in our version of libgcc. A bug has been filed against libgcc.

Concerning performance testing, we ran our 14 implementations on each of the entries of our result test vectors, timing the function call time in terms of machine cycles using the 1 and rounding mode combinations, as IEEE754-2008 specifies separate current rounding modes for the binary and decimal operations' subsets [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF] Intel rdtsc instruction after serialization with cpuid. We deducted the measured calling-time for an empty function with the same prototype. We ran the timing measurements repeatedly and with preheated caches, so to smooth out effects of modern superscalar processors. We made sure that the test vectors contained a representative, random subset of floating-point numbers, i.e. we made sure that special cases, such as addition and multiplication of zeros, NaNs etc. were represented but not over-represented. We present our timing measurements under the form of histograms which relate a certain cycle count interval to the number of times this cycle count interval was encountered. We performed the same performance testing on our optimized FMA implementation and on our GMP-based mixed-radix FMA one. For the sake of succinctness, we report only the performance testing results of one FMA function, called DBBB. This function takes all its arguments as IEEE754 binary64 numbers and returns a decimal64 number.

On the left of Figure 3 is the cycle count histograms for our DBBB FMA functions. It is easy to see that our implementations handle most cases in less than 50 cycles. Our implementations present a maximum call time of about 375 cycles, corresponding to cases when the rounding recovery phase gets used. As our algorithm is statically bounded in terms of memory consumption and makes no dynamic memory allocation, these maximum timing counts are hard bounds.

On the right of Figure 3 is the cycle count histograms of the GMP-based implementations for DBBB. For most cases, the GMP-based implementations is about 10 times slower than the implementation based on our algorithm presented in this article. The GMP-based maximum call times reach up to 46000 cycles in very rare cases, which is about 100 times slower than ours. This is due to the cost incurred by GMP's dynamical memory allocation mechanism.

V. CONCLUSION

In this article, we have presented an algorithm to compute a correctly-rounded IEEE754 binary64 or decimal64 result to a mixed-radix Fused-Multiply-and-Add operation, taking any combination of IEEE754 binary64 and decimal64 numbers in argument. Our algorithm trivially allows for covering correctly rounded mixed-radix addition, subtraction and multiplication and is designed in a way that allows mixed-radix division and square root to be computed with that FMA, too.

Our algorithm is based on a mixture of exact significand multiplication, finely precision-tuned radix conversion, binary addition and an exact rounding recovery phase that is launched if rounding of an approximation would not allow a correctly rounded result to be computed. The recovery phase is based on the observation that all IEEE754 binary64 and decimal64 floating-point numbers and their mid-points are integer multiples of a fixed 2 Q • 5 W , with constant Q, W ∈ Z, Q, W ≤ 0.

Our implementation can handle most input cases in about 50 machine cycles, which is a reasonable performance timing for a software-implemented floating-point operation. When the recovery phase gets used, latency stays below 400 cycles. Our implementation does not require any dynamic memory allocation, hence being suitable for integration into environments where dynamic memory management is not available.

Our implementation has been thoroughly tested against two other reference implementations, which we developed. Besides an issue with an undue IEEE754 inexact flag setting, which is outside our control, it tests correct. A detailed pencil-andpaper proof for our algorithm has been written and will be made available under the form of a research report.

While we are convinced that our implementation is one of the first step to paving the road to generalized mixed-radix floating-point arithmetic in IEEE754, quite a few points need to be addressed by future work: we subsumed all 14 possible mixed-radix binary and decimal input/output format combinations for the FMA under one single core algorithm. This eased development and the correctness proof; but it induces certain overhead due to exponent range over-estimation. Future work might perform a more detailed one-by-one analysis.

We showcased our mixed-radix FMA for the IEEE754 formats binary64 and decimal64 only. All possible combinations taking into account also the binary32 and decimal128 formats need to be covered by future development.

The necessity to fall back to exact arithmetic with pretty large integer-based accumulators for hard-to-round cases is not very intellectually satisfactory. Future work should try to address the issue of precomputing worst-case precision bounds for mixed-radix FMA again.

  double a=2000.0, c; _Decimal64 b=0.001D; c = a * b;

3 /

 3 * d = a * b + c * / 4 d = __builtin_fma(a,(double)b,c);

Fig. 1 .

 1 Fig. 1. Correctly Rounded Mixed-Radix FMA Algorithm

Fig. 2 .

 2 Fig. 2. Far-Path Addition/Subtraction Algorithm

Lemma 1 .- 2 X • 5 Y

 125 with 2 54 ≤ |t| < 2 55 and 2 109 ≤ |s| < 2 110 . Let us define X and Y such as X = N c -L and Y = P c -M ; X, Y ∈ Z. With this definition and the bounds of M , L, P c and N c given in equations (3) and (4), we can find the bounds of X and Y . For all X, Y ∈ Z with -3068 ≤ X ≤ 3230 and -1191 ≤ Y ≤ 1227, for s, t ∈ Z with with 2 54 ≤ |t| < 2 55 and 2 109 ≤ |s| < 2 110 s t ≥ η = 2 -177.61 .

  Finally we have |ε φ N P | ≤ 2 -63 . D. Rounding Test and recovery phase 1) Conversion to the output format: Once the multiplication and addition or subtraction are performed, we are left with a 64-bit binary significand with a binary exponent bounded such as 2 C • g with 2 63 ≤ |g| ≤ 2 64 -1 and -4174 ≤ C ≤ 3836.

  with ε * D the global error combining all the previous errors such as |ε * D | ≤ 2 -58.18 and 10 15 • 2 10 ≤ |q| ≤ (10 16 -1) • 2 10 . The rounding test boils down to compute the "magic bound" β D = ((10 16 -1) • 2 10 •

Fig. 3 .

 3 Fig. 3. Histogram of timing measurements for our implementation (left) and a reference implementation (right)

  represented by a table of 64-bit unsigned integers. The relative error arising during the computation of the table can be bounded as |ε τ0(µ) | ≤ 2 -63 .

  global error combining all the previous errors such as |ε * B | ≤ 2 -58.74 and 2 63 ≤ |g| ≤ 2 64 -1. The rounding test boils down to compute the "magic bound" β B = (2 64 • -10 • g 2 10 | ≥ β B . If this condition is not met, rounding to a 53-bit precision binary64 number is not possible. Therefore we compute f , the binary rounding boundary closest to 2

	ε * B B 1+ε *	and test if |g -2

  • 2 C • g into a normal binary64 value 2 A , another normal binary64 value 2 B • r and a normal or subnormal value 2 C • g such that a binary64 FMA computes 2 A ×2 B •r+2 C •g and rounds to the same binary64 value as (-1) s