Frédéric Goualard

𝑥/3) are not equal to 𝑥

Keywords: CCS Concepts: Mathematics of computing→Interval arithmetic, Theory of computation→Rounding techniques floating-point arithmetic, error analysis, roundoff, interval arithmetic iii

Rounding Error Analysis is routinely used to compute a worst-case error bound on the result of algorithms that use floating-point arithmetic. However, for some applications (e.g., when it is necessary to prove some inclusion of the result in a domain), the knowledge of both an upper-bound of the magnitude of the error and of its sign is paramount. Using standard rounding error analysis together with a simple systematic approach, we compute such information for the expressions 3𝑥, 𝑥/3, 3(𝑥/3), and 3𝑥/3, which can be used, e.g., in the proof of interval arithmetic operators.

Introduction

The implementation of Interval Arithmetic operators requires the evaluation of floating-point expressions, which can be marred by rounding errors that jeopardize the very properties that Interval Arithmetic is supposed to ensure [START_REF] Goualard | How do you compute the midpoint of an interval?[END_REF]. Rounding error analysis [START_REF] Nicholas | Accuracy and Stability of Numerical Algorithms[END_REF] can be used to determine in a systematic way the magnitude of these errors. It is usually, however, not concerned with their sign, a fatal flaw whenever the question of the inclusion of a value into an interval is to be ascertained. Some clever use of the binary floating-point number format properties as guaranteed by the IEEE 754 standard [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF] allows in special cases to work out the conditions under which rounding errors cancel out in such a way that the computed result is exactly equal to the true result. Several examples of such an approach can be found in Goldberg's survey [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF] and Muller et al.'s book [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], to name a few sources. On the other hand, few works concern themselves with determining both the magnitude and the sign of rounding errors in computing arithmetic expressions.

In this paper, we investigate the case of the expressions 3𝑥/3 and 3(𝑥/3) for 𝑥 a binary floating-point number. Using a mixture of classical Rounding Analysis and systematic study of the elementary operations "3𝑥" and "𝑥/3", we determinate for each possible value of 𝑥 the sign and the magnitude of the error in computing 3𝑥/3 and 3(𝑥/3). This study is used in a forthcoming article [START_REF] Goualard | Polyadic Splitting of FLoating-point Intervals[END_REF] to prove inclusion properties of some interval trisection operators. We expect it to be useful in the study of the many other algorithms that rely on a multiplication or a division by 3 (see, e.g., de Dinechin et al. 's work [START_REF] De | Table-Based Division by Small Integer Constants[END_REF]).

Binary floating-point numbers

In this study, we only concern ourselves with binary floating-point numbers as defined by the 1985 IEEE 754 standard [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF], by far the most widespread standard in use today for floating-point numbers. This paper assumes a prior knowledge of the basics of the IEEE 754 standard even though the most important points are summarized hereafter as needed for the sake of completeness.

We consider binary floating-point numbers from a set 𝔽 represented with a significand of size 𝑝 bits:

𝑥 ∈ 𝔽, 𝑥 = ±𝑏 0 .𝑏 1 … 𝑏 𝑝-1 × 2 𝑒
with an exponent 𝑒 in the range [𝐸 min , 𝐸 max ].

For the purpose of the systematic study of the elementary operations, we assume 𝑝 ⩾ 5, not a strong constraint considering that the smallest format presented in the IEEE 754 standard has 𝑝 = 24. We will not take into account the possibility of overflow when performing the multiplication by 3; on the other hand, the possibility of underflow will be fully adressed, albeit separately from the normal cases.

When the base 𝑏 used to represent a number 𝑥 is not unambiguously drawn from the context, we will indicate it as such: 𝑥 𝑏 .

All results are supposed to be rounded to nearest-even, the usual default rounding strategy. Given 𝑟 a real number, we note fl(𝑟) the floating-point value corresponding to 𝑟 rounded to nearest-even. Given 𝑥, a floating-point number, we note 𝑥 -the greatest floating-point number smaller than 𝑥, and 𝑥 + the smallest floating-point number greater than 𝑥. For an expression of the form (𝑥 ◇ 𝑦) □ 𝑧-with (𝑥, 𝑦, 𝑧) ∈ 𝔽 3 and (◇, □) arithmetic operators-, we use the shorthand fl⟨(𝑥 ◇ 𝑦) □ 𝑧⟩ 1 to mean fl ( fl ( 𝑥 ◇ 𝑦 ) □ 𝑧 ) . Lastly, we will extensively use in this paper the fact that rounding is order-preserving (rounding monotonicity [6, p. 38]):

∀(𝑟 1 , 𝑟 2 ) ∈ ℝ 2 ∶ 𝑟 1 ⩾ 𝑟 2 ⟹ fl(𝑟 1 ) ⩾ fl(𝑟 2 ) (1)
Note that if there is a strict inequality between 𝑟 1 and 𝑟 2 (i.e., 𝑟 1 > 𝑟 2 ), we still have fl(𝑟 1 ) ⩾ fl(𝑟 2 ).

With correctly rounded floating-point operators, as specified by the IEEE 754 standard, there is a simple relationship between a computed value and the real result [START_REF] Nicholas | Accuracy and Stability of Numerical Algorithms[END_REF]:

fl ( 𝑥 ◇ 𝑦 ) = (𝑥 ◇ 𝑦)(1 + 𝛿) + 𝜂, ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (𝑥, 𝑦) ∈ 𝔽 2 , 𝛿𝜂 = 0, |𝛿| ⩽ 𝑢, |𝜂| ⩽ 𝜇/2 ◇ ∈ {+, -, ×, ÷} (2)
1 Note the parentheses replaced by chevrons.

1

where 𝑢 is the unit roundoff equal to 2 -𝑝 and 𝜇 is the smallest positive subnormal number equal to 𝑢2 𝐸 min +1 . The bound on 𝛿 may be somewhat refined depending on the operator; in particular, we will make use of the following improved bounds [START_REF] Jeannerod | On relative errors of floating-point operations: optimal bounds and applications[END_REF]:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ fl ( 𝑥𝑦 ) = 𝑥𝑦(1 + 𝛿) , |𝛿| ⩽ 𝑢 1 + 𝑢 fl ( 𝑥 𝑦 ) = 𝑥 𝑦 (1 + 𝛿) , |𝛿| ⩽ 𝑢 -2𝑢 2 , (𝑥, 𝑦) ∈ 𝔽 2 , (3)
provided no underflow or overflow occurs.

Multiplying 𝑥 by three

Barring overflow, the multiplication of two binary floating-point numbers:

𝑥 = (-1) 𝑠 𝑥 𝑚 𝑥 × 2 𝑒 𝑥 and 𝑦 = (-1) 𝑠 𝑦 𝑚 𝑦 × 2 𝑒 𝑦
is given by [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]:

𝑥𝑦 = (-1) 𝑠 𝑥 ⊕𝑠 𝑦 𝑚 𝑥 𝑚 𝑦 × 2 𝑒 𝑥 +𝑒 𝑦 (4)
The multiplication of a floating-point number 𝑥 by 3 (i.e., 1.1 × 2 1 in normal floating-point binary form) is particularly easy to perform thanks to the simplicity of the representation of 3. We will consider first the case for which no underflow occurs; denormal numbers will be considered next.

The normal case

In this section, 𝑥 is a normal number of the form 𝑥 = (-1) 𝑠 𝑥 1.𝑏 1 … 𝑏 𝑝-1 × 2 𝑒 𝑥 . From Eq. (4), we see that we may consider 𝑥 > 0 only, as the results thus obtained need only be mirrored to obtain their counterpart when 𝑥 < 0.

Depending on 𝑥, there are two possibilities: either 3𝑥 has a 𝑝 + 1 bits significand (Eq. (5a)), or it has a 𝑝 + 2 bits significand (Eq. (5b)).

1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 ×2 𝑒 𝑥 × 1. 1 ×2 1 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 1. 𝑏 ′ 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 ×2 𝑒 𝑥 +1 (5a) 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 ×2 𝑒 𝑥 × 1. 1 ×2 1 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 1. 𝑏 ′ 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 𝑏 ′ 𝑝+1 ×2 𝑒 𝑥 +2 (5b)
From Eqs. (5a) and (5b), we can deduce conditions on 𝑥 for 𝑦 = 3𝑥 to have a 𝑝 + 1 bits significand; similarly, we can deduce the value of 𝑏 ′ 1 depending on the number of bits of 3𝑥:

Lemma 1. Given 𝑥 = 1.𝑏 1 … 𝑏 𝑝-1 × 2 𝑒 𝑥 a normal binary floating-point number. If 𝑦 = 3𝑥 has a 𝑝 + 1 bits significand, then 𝑥 = 1.0𝑏 2 … 𝑏 𝑝-1 × 2 𝑒 𝑥 and 𝑦 = 1.1𝑏 ′ 2 … 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 × 2 𝑒 𝑥 +1
Proof. Consider Eq. (5a) with the leftmost carries 𝑐 0 and 𝑐 1 added:

1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 ×2 𝑒 𝑥 × 1. 1 ×2 1 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 𝑐 0 𝑐 1 … 1. 𝑏 ′ 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 ×2 𝑒 𝑥 +1
We have:

1 + 𝑐 0 = 1 1 + 𝑏 1 + 𝑐 1 = 𝑐 0 𝑏 ′ 1 (henceforth, "𝑏 𝑖 𝑏 𝑗 " is to be interpreted as 𝑏 𝑖 × 2 + 𝑏 𝑗 )
Hence, 𝑐 0 = 0 and then

1 + 𝑏 1 + 𝑐 1 = 𝑏 ′ 1 , which can hold only if 𝑏 1 = 𝑐 1 = 0 and 𝑏 ′ 1 = 1. ⊓ ⊔ Lemma 2. Given 𝑥 = 1.𝑏 1 … 𝑏 𝑝-1 × 2 𝑒 𝑥 a normal binary floating-point number. If 𝑦 = 3𝑥 has a 𝑝 + 2 bits significand, then 𝑦 = 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 𝑏 ′ 𝑝+1 × 2 𝑒 𝑥 +2 .
Proof. Consider Eq. (5b) with the leftmost carries added:

1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 ×2 𝑒 𝑥 × 1. 1 ×2 1 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 1. 𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 + 𝑐 0 𝑐 1 1. 𝑏 ′ 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 𝑏 ′ 𝑝+1 ×2 𝑒 𝑥 +2
We have 1 + 𝑐 0 = 1𝑏 ′ 1 , from which we deduce 𝑐 0 = 1 and 𝑏 ′ 1 = 0. ⊓ ⊔

When 𝑦 = 3𝑥 has a 𝑝 + 2 bits significand, we have to discard the last two bits to be able to represent the result in 𝑝 bits format. If these two bits are "10"-that is, the real result is halfway between two representable numbers-, we need to know the antepenultimate bit of 𝑦's significand to decide the sign of the error. Consequently, we must consider the three leftmost bits of 𝑥. Appendix A presents all possible cases for the multiplication by 3.

For each possible value of 𝑥 with respect to the five bits of its significand considered above (The implicit '1' bit before the radix point, the first bit to the right of the radix point, and the three rightmost bits), Tables 1 and2 present all possible values for fl(3𝑥). For homogeneity-required later on when we will investigate the expressions 3𝑥/3 and 3(𝑥/3)-, fl(3𝑥) is also displayed with the same five bits, the exponents being left implicit. The third column "Error" gives the sign of the error:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ EQ ∶ fl(3𝑥) = 3𝑥 GT ∶ fl(3𝑥) > 3𝑥 LT ∶ fl(3𝑥) < 3𝑥

The denormal case

Consider now the case in which 𝑥 = 0.𝑏 1 … 𝑏 𝑝-1 × 2 𝐸 min is a denormal number. There are two cases depending on the value of 𝑏 1 . For 𝑏 1 = 0, the following result holds: Lemma 3. Given a denormal floating-point number of the form 𝑥 = 0.0𝑏 2 … 𝑏 𝑝-1 × 2 𝐸 min , we have:

fl(3𝑥) = 3𝑥 (6)
3 of 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝑒 𝑥 x fl(3x) Error 1.1𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 EQ 1.1𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 LT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 EQ 1.1𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 LT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 111 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 LT Proof. Trivial. Consider: 0. 0 𝑏 2 ⋯ 𝑏 𝑝-1 ×2 𝐸 min × 1. 1 ×2 1 0. 0 𝑏 2 𝑏 3 ⋯ 𝑏 𝑝-1 + 0. 0 𝑏 2 𝑏 3 ⋯ 𝑏 𝑝-1 0 𝑏 ′ 1 . 𝑏 ′ 2 𝑏 ′ 3 ⋯ 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 ×2 𝐸 min (7) ⊓ ⊔
On the other hand, if 𝑏 1 ≠ 0, we have:

0. 1 𝑏 2 ⋯ 𝑏 𝑝-1 ×2 𝐸 min × 1. 1 ×2 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-1 + 0. 1 𝑏 2 ⋯ 𝑏 𝑝-1 𝑐 0 𝑐 1 𝑏 ′ 0 . 𝑏 ′ 1 𝑏 ′ 2 𝑏 ′ 3 ⋯ 𝑏 ′ 𝑝-1 𝑏 ′ 𝑝 ×2 𝐸 min +1
where 𝑐 0 and 𝑐 1 are, as before, the leftmost carries.

It is easy to see that 𝑐 0 and 𝑐 1 must be equal. Furthermore:

{ if 𝑐 0 = 𝑐 1 = 0, then 𝑏 ′ 0 = 0 and 𝑏 ′ 1 = 1 if 𝑐 0 = 𝑐 1 = 1, then 𝑏 ′ 0 = 1 and 𝑏 ′ 1 = 0
Multiplying 𝑥 by 3 can then lead to a 𝑝 bits significand (when 𝑏 ′ 0 = 0) or to a 𝑝 + 1 bits significand (when 𝑏 ′ 0 = 1). Appendix B displays the different results for all possible values of the rightmost bits2 . From Appendix B and Eq. ( 7), we can deduce the value of fl(3𝑥) obtained from rounding 3𝑥 to nearest-even (see Table 3).

Dividing 𝑥 by three

Similarly to the multiplication, we will first consider the case in which 𝑥 is a normal number and no underflow occurs; we will then investigate what happens when 𝑥 is a denormal number or when the result underflows.

For 𝑥 = (-1) 𝑠 𝑥 𝑚 𝑥 × 2 𝑒 𝑥 and 𝑦 = (-1) 𝑠 𝑦 𝑚 𝑦 × 2 𝑒 𝑦 , the quotient 𝑥/𝑦 is given by [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]:

𝑥/𝑦 = (-1) 𝑠 𝑥 ⊕𝑠 𝑦 𝑚 𝑥 /𝑚 𝑦 × 2 𝑒 𝑥 -𝑒 𝑦 (8)
To determine the sign of the error of fl ( 𝑥/𝑦 ) relative to 𝑥/𝑦, we may only consider the significands 𝑚 𝑥 and 𝑚 𝑦 . Consequently, we will only display the exponents when rendered necessary by the context. We also restrict ourselves to positive values for 𝑥, the negative case being easily deduced.

The normal case

As for the multiplication, we need only consider the first bit after the radix point and the last three bits of the significand to determine the result of the division by 3. The following three lemmas tell us how many bits to compute in order to obtain a 𝑝 bits significand correctly rounded to nearest-even.

Lemma 4. If 𝑥 = 1.1𝑏 2 … 𝑏 𝑝-1 × 2 𝑒 𝑥
is a normal floating-point number, then we need to compute 𝑝 + 1 bits to get a correctly rounded 𝑝 bits significand for 𝑦 = fl(𝑥/3); in addition, we have:

fl ( 𝑦 ) = 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-1 × 2 𝑒 𝑥 -1 (9)
provided no underflow occurs.

Table 3: Result and sign of the error for the multiplication by

3 of 𝑥 = 0.𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝐸 min x fl(3x) Error 0.0𝑏 2 ⋯ 𝑏 𝑝-1 𝑏 ′ 0 .𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-1 EQ 0.1𝑏 2 ⋯ 𝑏 𝑝-4 000 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 00 EQ 0.1𝑏 2 ⋯ 𝑏 𝑝-4 001 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 10 GT 0.1𝑏 2 ⋯ 𝑏 𝑝-4 010 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 110 EQ 0.1𝑏 2 ⋯ 𝑏 𝑝-4 011 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 00 LT 0.1𝑏 2 ⋯ 𝑏 𝑝-4 100 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 10 EQ 0.1𝑏 2 ⋯ 𝑏 𝑝-4 101 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 00 GT 0.1𝑏 2 ⋯ 𝑏 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 01 EQ 0.1𝑏 2 ⋯ 𝑏 𝑝-4 111 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-3 10 LT Proof.
Consider the process of dividing 𝑥 by 3 at its beginning:

0 1 0 ⋯ 11 ) 1 1 𝑏 2 … 𝑏 𝑝-1 1 1 0 𝑏 2 ⋯
Due to the leading 0 of the quotient, we need to compute one bit more to get 𝑝 bits in the rounded significand. Equation ( 9) is obvious from the division process shown above. ⊓ ⊔ Lemma 5. If 𝑥 = 1.0𝑏 2 … 𝑏 𝑝-1 × 2 𝑒 𝑥 is a normal floating-point number, then we need to compute 𝑝 + 2 bits to get a 𝑝 bit significand for 𝑦 = 𝑥/3.

Proof. Consider the process of dividing 𝑥 by 3 at its beginning:

0 0 1 𝑏 ′ 1 ⋯ 11 ) 1 0 𝑏 2 𝑏 3 … 𝑏 𝑝-1 1 0 𝑏 2 𝑟 0 𝑟 1 𝑏 3 ⋯
Due to the two leading 0 of the quotient, we need to compute two bits more to get 𝑝 bits in the significand of the result. ⊓ ⊔ Lemma 6. Given 𝑥 = 1.𝑏 1 … 𝑏 𝑝-1 × 2 𝑒 𝑥 a normal floating-point number and 𝑦 = 𝑥/3. There are three possible forms for the significand of 𝑦:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1.𝑏 ′ 1 … 𝑏 ′ 𝑝-2 0 1.𝑏 ′ 1 … 𝑏 ′ 𝑝-2 01 1.𝑏 ′ 1 … 𝑏 ′ 𝑝-2 10
where 𝑠 denotes an infinite repetition of the binary string 𝑠.

Proof. In the process of dividing 𝑥 by 3, there are only three possible remainders at each step: 0 2 , 1 2 , or 10 2 . Consequently, when 𝑥 = 1.1𝑏 2 … 𝑏 𝑝-1 × 2 𝑒 𝑥 , we only have the following three situations (the case 𝑥 = 1.0𝑏 2 … 𝑏 𝑝-1 × 2 𝑒 𝑥 is analogous):

0 1 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-2 0 … 11 ) 1 1 𝑏 2 … 𝑏 𝑝-1 ⋯ ↓ 0 0 0 1 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-2 0 1 … 11 ) 1 1 𝑏 2 … 𝑏 𝑝-1 ⋯ ↓ ↓ 1 0 0 1 0 1 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-2 1 0 … 11 ) 1 1 𝑏 2 … 𝑏 𝑝-1 ⋯ ↓ 10 0 ↓ 1 0 (10) ⊓ ⊔
For each triplet of rightmost bits of the significand of 𝑥, there are three cases to consider depending on the remainder after having taken into account the bit 𝑏 𝑝-4 . All the cases are presented in Section C.

Tables 4 and5 summarize the results for the whole section.

The denormal case

There are two possibilities for an underflow to occur when dividing 𝑥 by 3: either 𝑥 is already a denormal number, or 𝑥 is normal but the rounded result of the division is a denormal number.

If 𝑥 is a denormal number, we have:

0. 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 × 2 𝐸 min 1.1 × 2 1 ) 0. 𝑏 1 … 𝑏 𝑝-4 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 × 2 𝐸 min ⋯ ⋯ … … (11)
If 𝑥 is a normal number of the form 𝑥 = 1.𝑏 1 … 𝑏 𝑝-1 × 2 𝐸 min , we have, either:

0. 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 × 2 𝐸 min 1.1 × 2 1 ) 1. 0 … 𝑏 𝑝-4 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 × 2 𝐸 min ⋯ ⋯ … … (12)
or:

0. 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 × 2 𝐸 min 1.1 × 2 1 ) 1. 1 … 𝑏 𝑝-4 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 × 2 𝐸 min ⋯ ⋯ … … (13)
In all three cases, the result of the division cannot be scaled to the right since the exponent is the smallest possible. As a consequence, we need only compute 𝑝 bits to obtain fl(𝑥/3

) = 0.𝑏 ′ 1 … 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 × 2 𝐸 min .
x fl(x/3) Error 1.0𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 2 ⋯ 𝑏 𝑝-4 111 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT Table 5:
Result and sign of the error for the division by 3 of 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝑒 𝑥 with no underflow.

x fl(x/3)

Error 1.1𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 LT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 LT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 EQ 1.1𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 LT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 GT 1.1𝑏 2 ⋯ 𝑏 𝑝-4 111 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 GT 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 EQ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 LT
One can easily check that, for 𝑥 = 0.𝑏 1 ⋯ 𝑏 𝑝-1 × 2 𝐸 min and 𝑥 = 1.0𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝐸 min , the bits 𝑏 ′ 𝑝-3 , 𝑏 ′ 𝑝-2 , 𝑏 ′ 𝑝-1 must be equal to, respectively, 𝑏 ′ 𝑝-5 , 𝑏 ′ 𝑝-4 , 𝑏 ′ 𝑝-3 as obtained when dividing a normal number of the form 1.0𝑏 … 𝑏 𝑝-1 × 2 𝑒 -with 𝑒 > 𝐸 min -by 3. When 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝐸 min , the bits 𝑏 ′ 𝑝-3 , 𝑏 ′ 𝑝-2 , 𝑏 ′ 𝑝-1 must be equal to, respectively, 𝑏 ′ 𝑝-4 , 𝑏 ′ 𝑝-3 , and 𝑏 ′ 𝑝-2 when dividing by 3 a normal number of the form 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝑒 , with 𝑒 > 𝐸 min .

For an underflow to occur when dividing a normal number 𝑥 by 3, 𝑥 must be smaller or equal to (3𝜆) -with (3𝜆) -= 1.01 ⋯ 111 × 2 𝐸 min +1 -, where 𝜆 = 1.0 × 2 𝐸 min is the smallest positive normal floating-point number. The case in which the exponent of 𝑥 is 𝐸 min has already been dealt with above. When 𝑥 is of the form 1.0𝑏 … 𝑏 𝑝-1 × 2 𝐸 min +1 , it is easy to see in Eq. (14):

0 0. 1 ⋯ 𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 𝑏 ′ 𝑝-1 ×2 𝐸 min 1.1 × 2 1 ) 1. 0 𝑏 2 … 𝑏 𝑝-4 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 ×2 𝐸 min+1 ⋯ ⋯ … … (14)
that we need to compute 𝑝 + 1 bits in order to get a 𝑝 bits significand (barring the rounding process), and that the last three bits of the result are obtained following the same rules as for the case in which 𝑥 is of the form 1.1𝑏 … 𝑏 𝑝 × 2 𝑒 with 𝑒 > 𝐸 min + 1.

Composing the multiplication and division

As seen in Section 2, the correctness of the floating-point multiplication and division are precisely defined by IEEE 754. From Eq. (2) and Eq. (3), we may compute the worst case error for the composition of a multiplication and a division. In theory, there are eight cases depending on the order of the operations and on the outcome of each operation (underflow or not). Tables 8 and9 summarize the error analysis when computing either fl⟨(3𝑥)/3⟩ or fl⟨3(𝑥/3)⟩. The crossed out cells in the tables correspond to cases that cannot occur: Consider, for example, Table 8; if there is an underflow when multiplying 𝑥 by 3, there will surely be an underflow when dividing 3𝑥 by 3.

Table 8: Worst case error analysis for fl⟨ 3𝑥 3 ⟩, with:

|𝜂 1 | ⩽ 𝜇/2, |𝜂 2 | ⩽ 𝜇/2, |𝛿 1 | ⩽ 𝑢/(1 + 𝑢), |𝛿 2 | ⩽ 𝑢 -2𝑢 2 . X X X X X X X X X X 𝑦 = 3𝑥 𝑧 = 𝑦/3 Underflow No underflow Underflow 3𝑥 + 𝜂 1 3 + 𝜂 2 X X X X X X X 3𝑥 + 𝜂 1 3 (1 + 𝛿 2 )
No underflow

3𝑥(1 + 𝛿 1 ) 3 + 𝜂 2 3𝑥 3 (1 + 𝛿 1 )(1 + 𝛿 2 )
Table 9: Worst case error analysis for fl⟨3 𝑥 3 ⟩, with:

|𝜂 1 | ⩽ 𝜇/2, |𝜂 2 | ⩽ 𝜇/2, |𝛿 1 | ⩽ 𝑢 -2𝑢 2 , |𝛿 2 | ⩽ 𝑢/(1 + 𝑢). X X X X X X X X X X 𝑦 = 𝑥/3 𝑧 = 3𝑦 Underflow No underflow Underflow 3 ( 𝑥 3 + 𝜂 1 ) + 𝜂 2 3 ( 𝑥 3 + 𝜂 1 ) (1 + 𝛿 2 )
No underflow X X X X X X X X 3 ( 𝑥

3 (1 + 𝛿 1 ) ) + 𝜂 2 3 𝑥 3 (1 + 𝛿 1 )(1 + 𝛿 2 )

Multiplication and division with no underflow

Considering first the case in which no underflow occurs, we get, irrespective of the order of the operations considered:

fl⟨ 3𝑥 3 ⟩ = fl⟨3 𝑥 3 ⟩ = 𝑥(1 + 𝛿 𝑎 )(1 + 𝛿 𝑏 )
with |𝛿 𝑎 | ⩽ 𝑢/(1 + 𝑢) and |𝛿 𝑏 | ⩽ 𝑢 -2𝑢 2 . Equivalently, we have:

fl⟨ 3𝑥 3 ⟩ = fl⟨3 𝑥 3 ⟩ = 𝑥(1 + 𝛿) ( 15
)
with 𝛿 = 𝛿 𝑎 + 𝛿 𝑏 + 𝛿 𝑎 𝛿 𝑏 , which means that |𝛿| ⩽ 2(𝑢 -2𝑢 3 )/(1 + 𝑢), that is |𝛿| < 2𝑢. The relative distance between two consecutive floating-point numbers is contained in the range [𝑢, 2𝑢]; Consequently, the rounded value of 3𝑥/3 or 3(𝑥/3) may not round to 𝑥 exactly, but we have the guarantee that it must be at worst one of the immediate floating-point neighbors of 𝑥 (See Fig. 1).

Multiplication and division with underflow

As can be seen in Table 8 and Table 9, the worst case errors are different for fl⟨(3𝑥)/3⟩ and fl⟨3(𝑥/3)⟩ whenever an underflow occurs. We must, therefore, consider separately the two computations.

The case of fl⟨(3𝑥)/3⟩

Underflow on multiplication and division. When both the multiplication and division underflow, we get:

fl⟨ 3𝑥 3 ⟩ = 3𝑥 + 𝜂 1 3 + 𝜂 2 = 𝑥 + ( 𝜂 1 3 + 𝜂 2 ) , with |𝜂 1 | ⩽ 𝜇 2 , |𝜂 2 | ⩽ 𝜇 2 Hence, fl⟨ 3𝑥 3 ⟩ = 𝑥 + 𝜂, with |𝜂| ⩽ 2𝜇 3 ( 16
)
For the multiplication to underflow, 𝑥 must be a subnormal number. As shown in Figure 2, the absolute distance between two consecutive subnormals is 𝜇. Consequently, we deduce from Eq. ( 16)

that fl⟨ 3𝑥 3 ⟩ = 3𝑥 3 = 𝑥
if both the multiplication and the division underflow. Underflow on division but not on multiplication. If the multiplication does not underflow but the division does, we get:

fl⟨ 3𝑥 3 ⟩ = 3𝑥(1 + 𝛿 1 ) 3 + 𝜂 2 , with |𝛿 1 | ⩽ 𝑢 1 + 𝑢 , |𝜂 2 | ⩽ 𝜇 2 .
It is easy to check that for 𝑥 = 𝜆, we have fl⟨(3𝑥)/3⟩ = 𝜆. By monotonicity of rounding, we get that for all 𝑦 greater than 𝑥, we have fl⟨(3𝑦)/3⟩ ⩾ 𝜆. Hence, 𝑥 must be a subnormal number for the division to underflow when computing (3𝑥)/3. Therefore, the absolute distance between 3𝑥 and fl(3𝑥) must be smaller or equal to 𝜇 (See Figure 2), which leads to:

fl⟨ 3𝑥 3 ⟩ = 3𝑥 + 𝜂 3 3 + 𝜂 2 , with |𝜂 3 | ⩽ 𝜇, |𝜂 2 | ⩽ 𝜇 2 = 𝑥 + 5 6 𝜇
As in the previous case, we deduce:

fl⟨ 3𝑥 3 ⟩ = 3𝑥 3 = 𝑥
if the multiplication does not underflow but the division does.

The case of fl⟨3(𝑥/3)⟩

Underflow on division and multiplication. If both the division and the multiplication underflow, we get:

fl⟨3 𝑥 3 ⟩ = 3 ( 𝑥 3 + 𝜂 1 ) + 𝜂 2 , with |𝜂 1 | ⩽ 𝜇/2, |𝜂 2 | ⩽ 𝜇/2 (17)
Using monotonicity of rounding, if fl⟨3(𝑥/3)⟩ underflows, we must have 𝑥 ⩽ 𝜆. Therefore, fl(𝑥/3) is of the form 0.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 ×2 𝐸 min . According to Table 3, the multiplication is then without error, that is: 𝜂 2 = 0. Equation (17) simplifies to:

fl⟨3 𝑥 3 ⟩ = 𝑥 + 3𝜂 1 , with |𝜂 1 | ⩽ 𝜇/2
Hence,

|fl⟨3 𝑥 3 ⟩ -𝑥| ⩽ 3𝜇 2
Since the absolute distance between subnormal numbers is 𝜇, we have:

fl⟨3 𝑥 3 ⟩ ∈ {𝑥 -, 𝑥, 𝑥 + }
Underflow on division but not on multiplication. If the division underflows but the multiplication does not, we get:

fl⟨3 𝑥 3 ⟩ = 3 ( 𝑥 3 + 𝜂 1 ) (1 + 𝛿 2 ), with |𝜂 1 | ⩽ 𝜇/2, |𝛿 2 | ⩽ 𝑢/(1 + 𝑢) (18)
It is trivial to check that we must have 𝑥 ∈ [𝜆 + , (3𝜆) -], using monotonicity of rounding. For convenience, we will investigate the actual worst case error on three disjoint domains for 𝑥:

• If 𝑥 ∈ [𝜆 + , ( 3 2 𝜆) -] (i.e., 𝑥 ∈ [1.00 ⋯ 001 × 2 𝐸 min , 1.011 ⋯ 111 × 2 𝐸 min ]).
According to Table 6, we have:

𝑥 ∈ [𝜆 + , ( 3 2 𝜆)
- ] ⟹ fl ( 𝑥 3 ) = 0.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 × 2 𝐸 min
According to Table 3, we therefore have:

fl ( 3fl ( 𝑥 3 )) = 3fl ( 𝑥 3 ) which means: fl⟨3 𝑥 3 ⟩ = 3( 𝑥 3 + 𝜂 1 ), |𝜂 1 | ⩽ 𝜇/2
Therefore:

|fl⟨3 𝑥 3 ⟩ -𝑥| ⩽ 3𝜇 2
Since, the absolute distance between floating-point numbers is 𝜇 in the domain [𝜆 + , ( 3 2 𝜆) -], we have:

∀𝑥 ∈ [𝜆 + , (3𝜆/2) -]∶ fl⟨3 𝑥 3 ⟩ ∈ {𝑥 -, 𝑥, 𝑥 + } (19) • If 𝑥 ∈ [3𝜆/2, (2𝜆) -] (i.e. 𝑥 ∈ [1.10 ⋯ 000 × 2 𝐸 min , 1.11 ⋯ 111 × 2 𝐸 min ])
It is easy to check, using monotonicity of rounding, that:

𝑥 ∈ [3𝜆/2, (2𝜆) -] ⟹ fl⟨3 𝑥 3 ⟩ ∈ [3𝜆/2, ((2𝜆) -) -]
From Table 3, we get for 𝑦 a subnormal number:

fl ( 3𝑦 ) = 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-1 × 2 𝐸 min ⟹ fl ( 3𝑦 ) = 3𝑦 Hence: ∀𝑥 ∈ [3𝜆/2, (2𝜆) -]∶ fl⟨3 𝑥 3 ⟩ = 3( 𝑥 3 + 𝜂 1 ) = 𝑥 + 3𝜂 1 , |𝜂 1 | ⩽ 𝜇/2
From which we get:

|fl⟨3 𝑥 3 ⟩ -𝑥| ⩽ 3𝜇 2 Consequently, we have: ∀𝑥 ∈ [3𝜆/2, (2𝜆) -]∶ fl⟨3 𝑥 3 ⟩ ∈ {𝑥 -, 𝑥, 𝑥 + } (20) • If 𝑥 ∈ [2𝜆, (3𝜆) -] (i.e., 𝑥 ∈ [1.0 × 2 𝐸 min +1 , 1.01 ⋯ 111 × 2 𝐸 min +1 ]).
It is trivial to show, using monotonicity of rounding, that:

𝑥 ∈ [2𝜆, (3𝜆) -] ⟹ fl⟨3 𝑥 3 ⟩ ∈ [2𝜆, ((3𝜆) -) -]
In the interval [2𝜆, ((3𝜆) -) -], the absolute distance between two consecutive floating-point numbers is 2𝜇 (See Figure 2), which means that:

∀𝑦 ∶ fl ( 3𝑦 ) ∈ [2𝜆, ((3𝜆) -) -] ⟹ |fl ( 3𝑦 ) -3𝑦| ⩽ 𝜂 3 , |𝜂 3 | ⩽ 𝜇
Therefore, for 𝑥 ∈ [2𝜆, (3𝜆) -]:

fl⟨3 𝑥 3 ⟩ = 3( 𝑥 3 + 𝜂 1 ) + 𝜂 3 , |𝜂 1 | ⩽ 𝜇/2, |𝜂 3 | ⩽ 𝜇 = 𝑥 + 𝜂, |𝜂| ⩽ 5𝜇/2
Then:

|fl⟨3 𝑥 3 ⟩ -𝑥| ⩽ 5𝜇 2
The distance between consecutive floating-point numbers in the domain involved being 2𝜇, we eventually get:

∀𝑥 ∈ [2𝜆, (3𝜆) -]∶ fl⟨3 𝑥 3 ⟩ ∈ {𝑥 -, 𝑥, 𝑥 + } (21)
From Equations (19), (20), (21), we deduce that for 𝑥 ∈ [𝜆 + , (3𝜆) -], fl⟨3(𝑥/3)⟩ is equal to 𝑥 itself or one of its immediate predecessor or successor:

∀𝑥 ∈ [𝜆 + , (3𝜆) -]∶ fl⟨3 𝑥 3 ⟩ ∈ {𝑥 -, 𝑥, 𝑥 + } (22)

Putting it all together

Knowing the worst case error for each operation, and armed with the tables from Sections 3 and 4, we may now determine the possible signs of the error when computing fl⟨(3𝑥)/3⟩ and fl⟨3(𝑥/3)⟩. Consider for example the floating-point number 𝑥 of the form 1.1𝑏 2 … 𝑏 𝑝-4 001 × 2 𝑒 𝑥 and the expression (3𝑥)/3. Using Table 2, we get all the possible forms for the significand of 𝑦 = fl(3𝑥). They are-in the absence of underflow: For each of these four forms, we get from Table 4 all the possible forms for the significand of fl ( 𝑦/3 ) . For example, for the first form 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-4 001, we get: We can then discard from Eq. (24) all forms but the second one: 1.1𝑏 2 … 𝑏 𝑝-4 001. If we do the same work for all forms in Eq. (23), we will find that the only possible significand for fl⟨(3𝑥)/3⟩ is 1.1𝑏 2 … 𝑏 𝑝-4 001, that is the one of 𝑥 itself. Doing the work exemplified above for all possible significands for 𝑥 that do not lead to an underflow in either the division or the multiplication, we get Table 10. As proved in Section 5.2.1, should an underflow arise during the computation of fl⟨(3𝑥)/3⟩, it is not necessary to compute such a table since we then always have fl⟨(3𝑥)/3⟩ = 𝑥.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-4 001 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 … 𝑏 ′ 𝑝-4 111 (23)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1.0𝑏 ″ 2 … 𝑏 ″ 𝑝-
Computing the same table for the expression fl⟨3(𝑥/3)⟩ results in a table with 96-that is, 2 5 × 3-entries: given 𝑥 = 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-3 𝑏 𝑝-2 𝑏 𝑝-1 × 2 𝑒 𝑥 , fl⟨3(𝑥/3)⟩ may always be 𝑥 or one of its two immediate neighbours, whatever the values of 𝑏 0 , 𝑏 1 , 𝑏 𝑝-3 , 𝑏 𝑝-2 , and 𝑏 𝑝-1 . Therefore, it is useless to display the resulting table as it is uninformative.

Conclusion

Using both traditional error analysis and the tables computed in Sections 3 and 4, we have been able to compute refined error bounds in Section 5 for fl⟨3𝑥/3⟩ and fl⟨3(𝑥/3)⟩. At first sight (See Tables 11 and12), it seems that the expression 3𝑥/3 gives overall more accurate results than the expression 3(𝑥/3) whenever an underflow occurs, and that no expression is better than the other when no underflow occurs.

⩽ 𝜇, |𝛿| < 2𝑢. X X X X X X X X X X 𝑦 = 𝑥/3 𝑧 = 3𝑦 Underflow No underflow Underflow {𝑥 -, 𝑥, 𝑥 + } {𝑥 -, 𝑥, 𝑥 + } No underflow - 𝑥(1 + 𝛿)
This is only part of the story, however. Our detailed analysis reveals that only when the last three bits of the significand of 𝑥 are "010" or "110" may the expression 3𝑥/3 be rounded to a value different from 𝑥 (see Table 10), while the same situation may occur whatever the value of 𝑥's significand with the expression 3(𝑥/3) (refer to the end of the previous section). The expression 3𝑥/3 is therefore to be prefered to 3(𝑥/3) in all cases.

Tables 1, 2, 3, 4, 5, 6, and 7 are all interesting and useful in their own right whenever one needs to obtain more information on the division and multiplication by three of a floating-point number than merely a worst-case error bound. We expect to reuse the method used here to prove properties regarding the expressions 3𝑥/3 and 3(𝑥/3), which combines traditional error analysis with systematic study of the behaviour of equivalence classes of floating-point numbers, in proving the correctness of polyadic splitting operators for interval analysis [START_REF] Goualard | Polyadic Splitting of FLoating-point Intervals[END_REF].

The question of multiplying and dividing a number by another in floating-point arithmetic has been addressed several times already in slightly different settings: Kahan [3, thm. 7] showed that fl⟨(𝑚/𝑛) × 𝑛⟩ = 𝑚 for integers 𝑚 and 𝑛 such that |𝑚| < 2 𝑝-1 and 𝑛 is the sum of two powers of 2. Edelman [START_REF] Edelman | When is 𝑥 × (1/𝑥) ≠ 1?[END_REF] investigated the problem of finding the smallest positive floating-point 𝑥 such that 𝑥(1/𝑥) ≠ 1.

Obviously, the next step for us is to address the problem presented in this paper for other integer values than 3. That might prove indeed necessary for our work on polyadic splitting operators. Once the preliminary worst case error analysis of an expression is performed manually, a large part of the proofs can be performed programmatically, which departs from the more involved-though more powerful and flexible-approaches used by, e.g., Kahan and Edelman. Generalizing the approach to other similar problems of interest is a direction of future research.

A Multiplication by 3 involving no denormal number

The two columns below present the results of the multiplication by 3 of a positive floating-point number 𝑥 depending on the last three bits of its significand. Thanks to Lemmas 1 and 2, we know that we also need to consider the first bit after the radix point to know how many bits from the result to discard. The left column (resp. right column) corresponds to the cases for which we need to discard the last bit (resp. the last two bits) underlined. For convenience, we display only the significands and leave the exponents implicit.

⋯

𝑏 ′ 𝑝-3 𝑏 ′ 𝑝-2 0 1

B Multiplication by 3 involving denormal numbers

List of the two different results for all possible values of the rightmost bits, on the left for a 𝑝 bits significand, and on the right for a 𝑝 + 1 bits significand, where the bit to discard is underlined.

0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 0 × 1. 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 0 + 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 0 0. 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 𝑏 ′ 𝑝-3 0 0 0 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 × 1. 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 + 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 1. 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-2 0 0 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 1 × 1. 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 1 + 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 0 1 0. 1 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 𝑏 ′ 𝑝-3 0 1 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 1 × 1. 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 1 + 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 1 1. 0 𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-2 1 1 0. 1 𝑏 2 ⋯ 𝑏 𝑝-4 0 1 0 × 1. 1 0.

C Division by 3 involving no denormal number

The leftmost column considers 𝑥 = 1.0𝑏 2 ⋯ 𝑏 𝑝-1 ×2 𝑒 𝑥 and the rightmost column considers 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 ×2 𝑒 𝑥 .

For each triplet of possible last three bits of the significand of 𝑥, we consider the three possible remainders on 𝑏 𝑝-4 .

When 𝑥 ends in "000"

0

Figure 1 :

 1 Figure 1: Orange disks representing the possible values for fl⟨3𝑥/3⟩ or fl⟨3(𝑥/3)⟩ with rounding to nearest-even in the absence of underflow.

Figure 2 :

 2 Figure 2: Absolute distance between small positive floating-point numbers

Table 1 :

 1 Result and sign of the error for the multiplication by 3 of 𝑥 = 1.0𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝑒 𝑥

	x 1.0𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 2 ⋯ 𝑏 𝑝-4 111	fl(3x) 2 ⋯ 𝑏 ′ 1.1𝑏 ′ 𝑝-4 000 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111	Error EQ EQ EQ EQ EQ EQ GT GT GT GT GT GT EQ EQ GT GT GT GT LT LT LT LT LT LT EQ EQ EQ EQ EQ EQ GT GT GT GT GT GT EQ EQ LT LT LT LT LT LT LT LT LT LT

Table 2 :

 2 Result and sign of the error for the multiplication by

Table 4 :

 4 Result and sign of the error for the division by 3 of 𝑥 = 1.0𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝑒 𝑥 with no underflow

Table 6 :

 6 Result and sign of the error for the division by 3 of 𝑥 = 1.𝑏 1 ⋯ 𝑏 𝑝-1 × 2 𝐸 min or 𝑥 = 0.𝑏 1 … 𝑏 𝑝-1 × 2 𝐸 min with underflow. As shown in Eqs. (11), (12) and (13), 𝑏 ′ 1 = 1 if and only if 𝑥 = 1.1𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝐸 min

	x 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 000 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 001 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 010 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 011 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 100 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 101 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 110 𝑏 0 .𝑏 1 𝑏 2 ⋯ 𝑏 𝑝-4 111	fl(x/3) 1 ⋯ 𝑏 ′ 0.𝑏 ′ 𝑝-4 000 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 011 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 101 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 000 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 011 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 110 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 001 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 110 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 011 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 001 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 100 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 110 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 001 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 100 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 111 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 010 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 100 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 111 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 010 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 101 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 111 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 010 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 101 0.𝑏 ′ 1 ⋯ 𝑏 ′ 𝑝-4 000	Error EQ GT LT LT EQ GT GT EQ LT EQ GT LT LT EQ GT GT LT EQ EQ GT LT LT EQ GT

Tables

6 and 7

summarize all the previous results for the denormal case.

Table 7 :

 7 Result and sign of the error for the division by 3 of 𝑥 = 1.0𝑏 2 ⋯ 𝑏 𝑝-1 × 2 𝐸 min +1 with underflow

	x 1.0𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 2 ⋯ 𝑏 𝑝-4 111	fl(x/3) 0.1 ⋯ 𝑏 ′ 𝑝-4 000 0.1 ⋯ 𝑏 ′ 𝑝-4 101 0.1 ⋯ 𝑏 ′ 𝑝-4 011 0.1 ⋯ 𝑏 ′ 𝑝-4 001 0.1 ⋯ 𝑏 ′ 𝑝-4 110 0.1 ⋯ 𝑏 ′ 𝑝-4 011 0.1 ⋯ 𝑏 ′ 𝑝-4 001 0.1 ⋯ 𝑏 ′ 𝑝-4 100 0.1 ⋯ 𝑏 ′ 𝑝-4 111 0.1 ⋯ 𝑏 ′ 𝑝-4 010 0.1 ⋯ 𝑏 ′ 𝑝-4 111 0.1 ⋯ 𝑏 ′ 𝑝-4 101 0.1 ⋯ 𝑏 ′ 𝑝-4 011 0.1 ⋯ 𝑏 ′ 𝑝-4 000 0.1 ⋯ 𝑏 ′ 𝑝-4 101 0.1 ⋯ 𝑏 ′ 𝑝-4 011 0.1 ⋯ 𝑏 ′ 𝑝-4 001 0.1 ⋯ 𝑏 ′ 𝑝-4 110 0.1 ⋯ 𝑏 ′ 𝑝-4 100 0.1 ⋯ 𝑏 ′ 𝑝-4 001 0.1 ⋯ 𝑏 ′ 𝑝-4 111 0.1 ⋯ 𝑏 ′ 𝑝-4 101 0.1 ⋯ 𝑏 ′ 𝑝-4 010 0.1 ⋯ 𝑏 ′ 𝑝-4 111	Error EQ LT GT GT EQ LT LT EQ GT EQ LT GT GT EQ LT LT GT EQ EQ LT GT GT EQ LT

Table 10 :

 10 Result and sign of the error when computing fl⟨(3𝑥)/3⟩ without underflow.

	x 1.0𝑏 2 ⋯ 𝑏 𝑝-4 000 1.0𝑏 2 ⋯ 𝑏 𝑝-4 001 1.0𝑏 2 ⋯ 𝑏 𝑝-4 010 1.0𝑏 2 ⋯ 𝑏 𝑝-4 011 1.0𝑏 2 ⋯ 𝑏 𝑝-4 100 1.0𝑏 2 ⋯ 𝑏 𝑝-4 101 1.0𝑏 2 ⋯ 𝑏 𝑝-4 110 1.0𝑏 2 ⋯ 𝑏 𝑝-4 111 1.1𝑏 2 ⋯ 𝑏 𝑝-4 000 1.1𝑏 2 ⋯ 𝑏 𝑝-4 001 1.1𝑏 2 ⋯ 𝑏 𝑝-4 010 1.1𝑏 2 ⋯ 𝑏 𝑝-4 011 1.1𝑏 2 ⋯ 𝑏 𝑝-4 100 1.1𝑏 2 ⋯ 𝑏 𝑝-4 101 1.1𝑏 2 ⋯ 𝑏 𝑝-4 110 1.1𝑏 2 ⋯ 𝑏 𝑝-4 111	fl⟨(3x)/3⟩ 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 010 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 110 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.0𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 000 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 001 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 011 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 100 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 101 1.1𝑏 ′ 2 ⋯ 𝑏 ′ 𝑝-4 111	Error EQ EQ GT EQ EQ EQ EQ EQ LT EQ EQ EQ GT EQ EQ EQ LT EQ

Table 11 :

 11 Refined worst case error analysis for fl⟨

	3𝑥 3 𝑧 = 𝑦/3 Underflow No underflow ⟩, with |𝛿| < 2𝑢. X X X X X X X X X X 𝑦 = 3𝑥
	Underflow	𝑥		-
	No underflow	𝑥		𝑥(1 + 𝛿)
	Table 12: Refined worst case error analysis for fl⟨3	𝑥 3	⟩, with: |𝜂|

Note that, even though only the last two bits need to be considered, we take into account the last three as in the previous section.

When 𝑥 ends in "011"

When 𝑥 ends in "100"

D Division by 3 involving denormal numbers

.

As explained in Section 4.2, the rules for the division involving denormal numbers may be inferred from the ones for the division involving normal numbers (Refer to Section C above).

Research Report

Rapport de recherche

When (3𝑥)/3 and 3(𝑥/3) are not equal to 𝑥 RR n°xx.xx Frédéric Goualard