10.6092/issn.1972-5787/4359
Double rounding occurs when a floating-point value is first rounded to an intermediate precision before being rounded to a final precision. The result of two such consecutive roundings can differ from the result obtained when directly rounding to the final precision. Double rounding practically happens, for instance, when implementing the IEEE754 binary32 format with an arithmetic unit performing operations only in the larger binary64 format, such as done in the PowerPC or x87 floating-point units. It belongs to the folklore in the floating-point arithmetic community that double rounding is innocuous for the basic arithmetic operations (addition, division, multiplication, and square root) as soon as the final precision is about twice larger than the intermediate one. This paper adresses the formal proof of this fact considering underflow cases and its extension to radices other than two.
application/pdf
1972-5787
4359
1
en
Pierre Roux
Journal of Formalized Reasoning
Coq
Innocuous Double Rounding of Basic Arithmetic Operations
Text.Serial.Journal
https://jfr.unibo.it/
7
Articles
2014-11-14
2014-06-27
2014-07-17
2016-01-29
131-142
JFR
Pierre Roux
2014/11/14
10.6092/issn.1972-5787/4359
131
ISAE, ONERA
1972-5787
1
Coq
en
142
Innocuous Double Rounding of Basic Arithmetic Operations
https://jfr.unibo.it/article/download/4359/4017
7
Innocuous Double Rounding
of Basic Arithmetic Operations
Pierre Roux
ISAE, ONERA


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




1.   INTRODUCTION
Floating-point numbers are commonly used to efficiently perform numerical com-
putations, which are then performed with a bounded precision. That is, only a
finite number of bits are kept after each arithmetic operation. Multiple choices are
available for this precision, typical examples being 24 bits for the IEEE754 binary32
format and 53 for the binary64 format.
   Double rounding occurs when a value is first rounded to an intermediate preci-
sion before being rounded to a final precision. The result of two such consecutive
roundings can differ from the result obtained when directly rounding to the final
precision. In radix 10, for instance, 1.495 could be rounded to 1.50 then to 2 if first
rounding to a 3-digit precision then to one digit whereas a direct rounding to the
nearest value with one digit would give 1. Double rounding practically happens,
for instance, when implementing the IEEE754 binary32 [IEE08] format with an
arithmetic unit only performing operations in the larger binary64 format, such as
done in the PowerPC or x87 floating-point units.
   It belongs to the folklore in the floating-point arithmetic community that double
rounding is innocuous for the basic arithmetic operations (addition, subtraction,
multiplication, division and square root) as soon as the final precision is about
twice larger than the intermediate one. The IEEE754 binary32 and binary64 for-
mats [IEE08] fulfill this condition. Figueroa [Fig95] published pen-and-paper proofs
of these results but only considering radix 2 and completely ignoring underflows. We
try to give more general results in this paper. Moreover, such proofs tend to involve
quite a number of rather subtle corner cases which make them particularly error
prone. We therefore use Coq [BCHPM04, Coq12] to ensure no such corner case has
been missed. The impact of double rounding on some particular floating-point al-

This work was done while the author was a visiting researcher at LRI, Inria Saclay – Île-de-France.

                                      Journal of Formalized Reasoning Vol. 7, No. 1, 2014, Pages 131–142.
132     ·      Pierre Roux

gorithms such as the Veltkamp/Dekker’s algorithms has already been studied using
Coq [MDMM13], but its innocuousness, under some conditions, on basic arithmetic
operations is only briefly mentioned.
  More formally, given a floating-point format F, a rounding ◦1 in this format and
a second rounding ◦2 in a larger-precision format, we want to ensure:
                         ∀x, y ∈ F, ◦1(◦2(x  y)) = ◦1(x  y)
                                               √
under some conditions, for  ∈ {+, −, ×, /, }. That is, the second rounding ◦2
is innocuous for basic arithmetic operations. This paper aims at studying under
which conditions the previous holds and to formally prove such results using a proof
assistant (Coq). It is worth noting that the results are highly dependent on the
                                                                                 √
operation  and that they may not hold for operations other than {+, −, ×, /, }.
Since all the proofs presented in the remaining of the paper have been formally
verified within Coq, we only present proof sketches. Although not all subcases are
precisely demonstrated, enough is shown to give a good idea of what is going on.
   The case of directed roundings (toward −∞ or +∞ or roundings to zero or
away) is rather simple: the above result holds if and only if ◦1 and ◦2 are rounding
in the same direction. The remaining will therefore focus on the case where ◦1
and ◦2 are roundings to nearest (possibly with different tie-break rules). We can
nevertheless notice that the results still hold (under possibly weaker hypotheses) if
◦1 is a directed rounding and ◦2 a rounding to nearest. The converse case (◦1 to
nearest and ◦2 directed) however usually requires slightly stronger hypothesis.
   Finally, we do not deal with overflows since they are harmless as soon as any
number which can be represented in the smaller precision format F can also be
represented in the larger precision, which is a reasonable assumption on exponent
ranges.
   The remainder of this section introduces the floating-point formalisms used in our
proofs. Then, Sections 2, 3, 4, and 5 illustrate our proofs about respectively the
multiplication, the addition/subtraction, the square root and the division. For the
sake of clarity, the proofs are presented on a formalism not modeling underflows,
although our actual Coq proofs handle them and precise requirements on underflow
will be given. Eventually, Section 6 summarizes the results and concludes.
   All our Coq developments are part of the Flocq library [BM11] available1 at
http://flocq.gforge.inria.fr/.
1.1   Floating-point Formats
We first formally define some floating-point formats used throughout the remaining
of the paper. A generic format is first defined, then particular cases: a fixed
precision format without underflow, a format with gradual underflows and one
with abrupt underflows. All these formats are formally defined in the Coq library
for floating-point arithmetic Flocq [BM11].
  For this purpose, we need a discrete logarithm.

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

with version 2.4.0 of Flocq.




Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
                 Innocuous Double Rounding of Basic Arithmetic Operations          ·     133

  Definition 1. Given β ∈ Z, β ≥ 2, we define the discrete logarithm lnβ (x) for
radix β of a real number x ∈ R \ {0} as the unique value ex in Z satisfying
                                β ex −1 ≤ |x| < β ex .
  A generic format Fϕ is defined by a function ϕ : Z → Z giving, for all numbers
sharing the same discrete logarithm, the precision they must be encoded with.
  Definition 2. Given ϕ : Z → Z, a value x ∈ R is said to be in the generic format
Fϕ if x = 0 or if there is some m ∈ Z such that
                                 x = m β ϕ(lnβ (x)) .
m is then called the mantissa of x and ϕ(lnβ (x)) its canonical exponent. In the
Flocq library, the proposition generic format β ϕ x means that x ∈ Fϕ .
  This generic format then allows us to define more concrete formats. Given some
p ∈ Z, p ≥ 1, the format FLXp models precision p floating-point numbers without
underflow.
   Definition 3. Given p ∈ Z, p ≥ 1, the format FLXp is defined as Fϕ with ϕ :
e 7→ e − p (denoted FLX exp p in Flocq).
   Although this format remains mostly theoretical since hardware floating-point
numbers can underflow, it is still a good model of what happens with actual floating-
point values when no underflow occurs.
   Given p ∈ Z, p ≥ 1, and emin ∈ Z, the format FLTp,emin models precision p
floating-point numbers with underflow handled by subnormals of exponent emin .
  Definition 4. Given p, emin ∈ Z, p ≥ 1, the format FLTp,emin is defined as Fϕ
with ϕ : e 7→ max(e − p, emin ) (denoted FLT exp emin p in Flocq).
  Remark 5. FLT24,−149 with radix β = 2 corresponds to the IEEE754 binary32
format and FLT53,−1074 to binary64 (neglecting overflows).
   Finally, given p ∈ Z and emin ∈ Z, the format FTZp,emin models precision p
floating-point numbers with underflow handled by flushing either to zero or to the
smallest normal number.
  Definition 6. Given p, emin ∈ Z, the format FTZp,emin is defined as Fϕ with
ϕ : e 7→ e − p when e ≥ emin + p and e 7→ emin + p − 1 otherwise (denoted
FTZ exp emin p in Flocq).
  The ϕ functions defining FLX, FLT and FTZ formats are illustrated on Figure 1.
  It only remains to define roundings.
 Definition 7. A rounding on a floating point format F is a function ◦ : R → F
monotone and equal to the identity on F, that is
                          ∀x, y ∈ R, x ≤ y ⇒ ◦(x) ≤ ◦(y)
                          ∀f ∈ F, ◦(f ) = f.
This corresponds to the hypothesis Valid rnd rnd in Flocq while the rounding
function ◦ on Fϕ is denoted by round β ϕ rnd.

                                         Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
134      ·     Pierre Roux


ϕ(e)                               ϕ(e)                           ϕ(e)
                              p                              p                               p


                                                                 emin +
                                                                  p−1

                                  emin                           emin


                              e             emin emin + p   e              emin emin + p     e

         (a) FLXp                        (b) FLTp,emin                 (c) FTZp,emin

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


   Definition 8. A rounding to nearest is a rounding ◦ that returns a nearest value
in F:
                             ∀x ∈ R, ∀f ∈ F, |x − ◦(x)| ≤ |x − f | .
They are of the form round β ϕ (Znearest ) in Flocq.
  Remark 9. In case x lies halfway between two consecutive values in F, ◦(x) is
not uniquely defined as it can be any of these two consecutive values. This choice is
called a tie-break rule. Since most of our proofs are valid regardless of the tie-break
rule (i.e., for any rounding to nearest) we don’t extend on this point.

2.    MULTIPLICATION
Proofs for the addition/subtraction being rather involved, let us begin with the
multiplication for which proofs are much simpler.
2.1    General Case
Figueroa [Fig95] proved double rounding to be innocuous for the multiplication if
it is performed with a precision at least twice as large. Although his proof only
addressed radix 2, this holds for any radix.
   Theorem 10. ([Fig95], double round mult FLX in our Coq development) For
p1 , p2 ∈ Z, if p2 ≥ 2p1 , then for ◦1 and ◦2 any roundings (for instance, toward −∞
or to nearest with any tie) respectively in FLXp1 and FLXp2 :
                         ∀x, y ∈ FLXp1 , ◦1(◦2(x × y)) = ◦1(x × y) .
     This in particular applies to binary32/64 when no underflow occurs.
   Remark 11. [Fig95] This bound 2p1 is optimal as shown by the following coun-
terexample. With radix β = 2, for p1 = 4, p2 = 7, we have for x = y = 1.1012
             ◦1(◦2(x × y)) = ◦1(◦2(10.1010012 )) = ◦1(10.101002 ) = 10.102
which differs from ◦1(10.1010012 ) = 10.112 for ◦1 and ◦2 appropriate roundings to
nearest (for instance, the common rounding to nearest with tie-break to even).
Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
                   Innocuous Double Rounding of Basic Arithmetic Operations             ·     135

  Proof. x and y having p1 significand digits, x × y has at most 2p1 ≤ p2 digits
hence ◦2(x × y) = x × y.
  The previous result still applies to binary32/64 even in case of underflow. Indeed,
proofs for the FLX, FLT and FTZ formats (double round mult {FLX,FLT,FTZ} in
our Coq development) are just corollaries of a single proof (double round mult)
carried out on the generic format Fϕ . This general proof shows innocuity of double
rounding of multiplication under the assumption
 ∀ex , ey ∈ Z, ϕ2 (ex + ey ) ≤ ϕ1 (ex ) + ϕ1 (ey ) ∧ ϕ2 (ex + ey − 1) ≤ ϕ1 (ex ) + ϕ1 (ey ),
the main argument being the same than in the above proof.
2.2   Odd Radix
It is worth noting that a much better result can be obtained for odd radices β, since
the hypothesis p2 ≥ 2p1 in Theorem 10 can be replaced with the weaker p2 ≥ p1
when β is odd. This nice result is, at least currently2 , certainly perfectly useless.
We nevertheless chose to develop it for the sake of exhaustivity and considering it
did not implied a huge overhead on our Coq development.
  Theorem 12. (double round mult beta odd FLX) When β is odd, for p1 , p2 ∈
Z, if p2 ≥ p1 , then for ◦1 and ◦2 roundings to nearest, with any tie, respectively in
FLXp1 and FLXp2
                       ∀x, y ∈ FLXp1 , ◦1(◦2(x × y)) = ◦1(x × y) .
  To give an idea of the proof, we first need a few additional definitions and lemmas.
  Definition 13. ◦↓ϕ denotes rounding towards −∞, i.e., ◦↓ϕ(x) is the largest floating-
point value in Fϕ below x:
                          ◦↓ϕ : x ∈ R 7→ max {f ∈ Fϕ | f ≤ x}.
  Definition 14. ulpϕ(x) denotes the unit in last place of x in format Fϕ :
                                   ulpϕ(x) := β ϕ(lnβ (x)) .
  Definition 15. midpϕ(x) denotes the following value, called midpoint
                                                 ulpϕ(x)
                              midpϕ(x) := ◦↓ϕ(x) +
                                                    2
as it lies halfway between two consecutive floating-point values around x.
  Lemma 16. (double round lt mid further place) Given two floating-point for-
mats Fϕ1 and Fϕ2 , two roundings to nearest, with any tie, ◦1 and ◦2 in these for-
mats, for all x ∈ R, if x > 0, ϕ2 (lnβ (x)) < ϕ1 (lnβ (x)), ϕ1 (lnβ (x)) ≤ lnβ (x)
and
                                             ulpϕ2(x)
                           x < midpϕ1(x) −
                                                2
then
                                     ◦1(◦2(x)) = ◦1(x) .
2 There have been a couple of computers using a radix 3 arithmetic, e.g., Setun built in Moscow
in the 60s and Ternac developped in New York in the 70s.
                                              Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
136        ·       Pierre Roux

                                                                           ◦2(x)

                                                                             x
   ◦↓ϕ1 (x)                                                                                             midpϕ1 (x)

                                                       ulpϕ (x)
                                                           1
                                                          2                                  ulpϕ (x)
                                                                                                 2
                                                                                                2


Fig. 2: Illustration of Lemma 16. Since x is below midpϕ1 (x), ◦1(x) = ◦↓ϕ1 (x) and since x is below
midpϕ1 (x) − ulpϕ2 (x) /2 then ◦2(x) is also below midpϕ1 (x), implying ◦1(◦2(x)) = ◦↓ϕ1 (x).

                                                ulpϕ (x1 )
                                                    2
                                                    2


      x0                                 x1                  m                   x2                           x3
                            ulpϕ (x0 )
                                1
                                2


Fig. 3: Illustration of Lemma 17. In radix β = 3, with ϕ2 = ϕ1 − 1, x0 and x3 are consecutive
values in Fϕ1 while x0 , x1 , x2 and x3 are consecutive values in Fϕ2 . The midpoint m in Fϕ1 is
                                          ulpϕ (x0 )              ulpϕ (x1 )
                                              1                       2
also a midpoint in Fϕ2 : m = x0 +             2
                                                       = x1 +         2
                                                                             .   This property relies on β being odd.


  This is illustrated on Figure 2.
  Proof. Since ϕ2 (lnβ (x)) ≤ ϕ1 (lnβ (x)), we get ◦2(x)−◦↓ϕ1(x) = ◦2(x) − ◦↓ϕ1(x) .
                                                                                                            ulpϕ2(y)
Moreover, ◦2 being a rounding to nearest, we get for all y ∈ R, |◦2(y) − y| ≤                                  2
                                     ulpϕ2(x)                                                     ulpϕ1(x)−ulpϕ2(x)
and the hypothesis x < midpϕ1(x)−       2     gives us x − ◦↓ϕ1(x) <                                      2         .
                                                      ulpϕ1(x)
Thus ◦2(x) − ◦↓ϕ1(x) ≤ |◦2(x) − x| + x − ◦↓ϕ1(x) <        2    . ◦1 being                           a rounding to
                                                        ↓
nearest, ◦1(◦2(x)) and ◦1(x) are then both equal to ◦ϕ1(x).
  Contrary to the previous lemmas which were valid for any radix β, the following
result shows a particularly interesting property of odd radices.
  Lemma 17. (midpoint beta odd remains) When β is odd, for all ϕ1 , ϕ2 : Z →
Z, for all x ∈ Fϕ1 , if ϕ2 (lnβ (x)) ≤ ϕ1 (lnβ (x)), then there exists an x0 ∈ Fϕ2 such
               ulpϕ1(x)            ulpϕ2(x)
that x +          2       = x0 +      2     .
  This result is illustrated on Figure 3. Intuitively, this means that, with odd
                               ulpϕ1(x)
radices, if a real number x +     2     is a midpoint for some precision, it is also
a midpoint in any larger precision. This is proved by induction on ϕ1 (lnβ (x)) −
ϕ2 (lnβ (x)) ∈ N and makes it possible to prove the following lemma.
  Lemma 18. (neq midpoint beta odd) When β is odd, for all ϕ1 , ϕ2 : Z → Z,
for all x ∈ R, if ϕ2 (lnβ (x)) ≤ ϕ1 (lnβ (x)) and x 6= midpϕ1(x) then ◦1(◦2(x)) = ◦1(x)
for any rounding to nearest, with any tie, ◦1 and ◦2 respectively in Fϕ1 and Fϕ2 .
  This means that with odd radices, double rounding is innocuous for all real values
which are not midpoints.



Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
                     Innocuous Double Rounding of Basic Arithmetic Operations           ·      137

                                                                ulpϕ2(x)
     Proof. Lemma 16 concludes when x < midpϕ1(x)−                 2     .   Let us then assume
                                       ulpϕ2(x)
                         midpϕ1(x) −            ≤ x < midpϕ1(x) .
                                          2
                                                                                   ulpϕ2(x)
From Lemma 17, there exist an x0 ∈ Fϕ2 such that x0 = midpϕ1(x) −                     2       and
                                                  ulpϕ2(x)
                                  x0 ≤ x < x0 +
                                                     2
hence ◦2(x) = x0 . Since ◦↓ϕ1(x) ≤ x0 < midpϕ1(x), then ◦1(x0 ) = ◦↓ϕ1(x) = ◦1(x).
The proof is similar for x > midpϕ1(x).
     A last lemma is then needed to prove the theorem.
  Lemma 19. (float neq midpoint beta odd) When β is odd, for all x ∈ R,
ϕ : Z → Z, if there exist m, e ∈ Z such that x = mβ e , then x 6= midpϕ(x)
     Proof. If e ≥ ϕ(lnβ (x)), then both x and ◦↓ϕ(x) are multiples of β ϕ(lnβ (x))
           ulpϕ(x)
whereas       2      is not. Otherwise, Lemma 17 allows us to conclude by a similar
reasoning.
   Proof (Theorem 12). Since x × y = (mx × my ) β ex +ey , with mx , my and
ex , ey the mantissas and exponents of x and y, Lemmas 18 and 19 make it possible
to conclude.

3.    ADDITION
3.1    General Case
Figueroa [Fig95] proved double rounding to be innocuous for the addition if it is
performed with a precision at least strictly twice as large. Although his proof only
addressed radix 2, this holds for any radix.
  Theorem 20. [Fig95] (double round plus FLX) For p1 , p2 ∈ Z, if p2 ≥ 2p1 +1,
then for ◦1 and ◦2 roundings to nearest, with any tie, respectively in FLXp1 and
FLXp2 :
                         ∀x, y ∈ FLXp1 , ◦1(◦2(x + y)) = ◦1(x + y) .
     This applies to IEEE754 binary32/64 when no underflow occurs.
  Remark 21. [Fig95] This bound 2p1 + 1 is optimal as shown by the following
counterexample. With radix β = 2, for p1 = 4, p2 = 8, we have for x = 1.0012 and
y = 0.000011112
             ◦1(◦2(x + y)) = ◦1(◦2(1.001011112 )) = ◦1(1.00110002 ) = 1.0102
which differs from ◦1(1.001011112 ) = 1.0012 when ◦1 and ◦2 are roundings to near-
est with tie-break to even.
     We need the following lemma to prove the previous theorem.
  Lemma 22. For all ϕ : Z → Z, x ∈ Fϕ , y ∈ R, if x > 0, y > 0 and lnβ (y) ≤
ϕ(lnβ (x)) − 2, then
                          0 < (x + y) − ◦↓ϕ(x + y) < β ϕ(lnβ (x))−2 .
                                              Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
138      ·     Pierre Roux

     Proof. Since lnβ (y) ≤ ϕ(lnβ (x)) − 2, we get y < ulpϕ(x) = ulpϕ(x + y), hence
◦↓ϕ(x  + y) = x, that is (x + y) − ◦↓ϕ(x + y) = y which implies the result.
  Proof (Theorem 20). We will only consider the case x > 0, y > 0 and y ≤ x,
other cases being pretty similar or easy to deduce. Either lnβ (y) ≥ ϕ1 (lnβ (x))−1 in
which case x+y has at most 2p1 +1 significand digits, which means that ◦2(x + y) =
x + y, or lnβ (y) ≤ ϕ1 (lnβ (x)) − 2 in which case the result follows from Lemmas 16
and 22.
  The previous result still applies to IEEE754 binary32/64 even in case of underflow
(double round plus FLT), with a similar proof thanks to the generic format Fϕ .
Indeed, proofs for formats with or without underflow are derived as immediate
corollaries of a proof on the generic format.
  Under the assumption β ≥ 3, the bound 2p1 + 1 can be replaced by 2p1 in
Theorem 20. The proof is similar (double round plus beta ge 3 FLX).
  Remark 23. This bound 2p1 is optimal as shown by the following counterex-
ample. With radix β = 10, for p1 = 4, p2 = 7, we have for x = 1.00110 and
y = 0.000499510
             ◦1(◦2(x + y)) = ◦1(◦2(1.001499510 )) = ◦1(1.00150010 ) = 1.00210
which differs from ◦1(1.001499510 ) = 1.00110 when ◦1 and ◦2 are roundings to
nearest with tie-break to even.
  Finally, we can notice that the same results can be immediately deduced for the
subtraction.
3.2     Odd Radix
As for the multiplication, a way better bound is obtained for odd radices.
  Theorem 24. (double round plus beta odd FLX) When β is odd, for p1 , p2 ∈
Z, if p2 ≥ p1 , then for ◦1 and ◦2 roundings to nearest, with any tie, respectively in
FLXp1 and FLXp2 :
                         ∀x, y ∈ FLXp1 , ◦1(◦2(x + y)) = ◦1(x + y) .
  Proof. Denoting mx , my and ex , ey the mantissas and exponents of x and y,
x+y = mβ e with e = min(ex , ey ) and m = mx β ex −min(ex ,ey ) +my β ey −min(ex ,ey ) ∈ Z.
Lemmas 18 and 19 make it possible to conclude.

4.    SQUARE ROOT
4.1     General Case
Figueroa [Fig95] proved double rounding to be innocuous for the square root if
it is performed with a precision larger than twice the original precision plus two.
Although his proof only addressed radix 2, this holds for any radix.
  Theorem 25. [Fig95] (double round sqrt FLX) For p1 , p2 ∈ Z, if p2 ≥ 2p1 +2,
then for ◦1 and ◦2 roundings to nearest respectively in FLXp1 and FLXp2
                                          √          √ 
                       ∀x ∈ FLXp1 , ◦1 ◦2 x = ◦1 x .
     This applies to IEEE754 binary32/64 when no underflow occurs.
Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
                  Innocuous Double Rounding of Basic Arithmetic Operations           ·     139

  Remark 26. [Fig95] This bound 2p1 + 2 is optimal as shown by the following
counterexample. With radix β = 2, for p1 = 4, p2 = 9, we have for x = 0.11112
            √ 
       ◦1 ◦2 x = ◦1(◦2(0.1111011111 . . .2 )) = ◦1(0.1111100002 ) = 1.0002
which differs from ◦1(0.1111011111 . . .2 ) = 0.11112 when ◦1 and ◦2 are roundings
to nearest with tie-break to even.
                                                          √
                               √              √      ulp ( x)
  Proof (Theorem 25). If x < midpϕ1( x)− ϕ22                  , the result follows from
                                                                                √
                                                 √               √      ulpϕ2( x)
Lemma 16. Similarly, the result holds when x > midpϕ1( x) +                   2    . We
will now show√  that  x cannot  lie
                                 √  between  those
                                                 √ two  bounds.    Let us   first denote
u1 := ulpϕ1( x), u2 := ulpϕ2( x), a := ◦↓ϕ1( x), b := u1 −u     2
                                                                   2
                                                                     and b0 := u1 +u2
                                                                                      2
                                                                                        .
Then assuming
                                         √
                                 a + b ≤ x ≤ a + b0
we get
                  a2 + u1 a − u2 a + b2 ≤ x ≤ a2 + u1 a + u2 a + b02 .
a2 +u1 a as well as x can be divided by u21 . However −u2 a+b2 > 0 and u2 a+b02 < u21
which is absurd.
  The previous result still applies to IEEE754 binary32/64 even in case of underflow
(double round sqrt FLT).
  Under the assumption β ≥ 4, the bound 2p1 + 2 can be replaced by 2p1 + 1 in
Theorem 25. The proof is similar (double round sqrt beta ge 4 FLX).
  Remark 27. This bound 2p1 + 1 is optimal as shown by the following counterex-
ample. With radix β = 10, for p1 = 4, p2 = 8, we have for x = 0.999910
           √ 
      ◦1 ◦2 x = ◦1(◦2(0.999949998 . . .10 )) = ◦1(0.9999500010 ) = 1.00010
which differs from ◦1(0.999949998 . . .10 ) = 0.999910 when ◦1 and ◦2 are roundings
to nearest with tie-break to even.
4.2   Odd Radix
As for the multiplication and addition, a way better bound is obtained for odd
radices.
  Theorem 28. (double round sqrt beta odd FLX) When β is odd, for p1 , p2 ∈
Z, if p2 ≥ p1 , then for ◦1 and ◦2 roundings to nearest respectively in FLXp1 and
FLXp2 :
                                          √        √ 
                        ∀x ∈ FLXp1 , ◦1 ◦2 x = ◦1 x .
                                                                √            √
  Proof. According to Lemma 18, it is enough to prove that x 6= midpϕ1( x).
Assuming the contrary, we get
                                            u2
                              x = r2 + ru +
                                            4
            ↓
               √                  √                                        √
where r := ◦ϕ1( x) and u := ulpϕ1( x). Then, when ϕ1 (lnβ (x)) ≥ 2ϕ1 (lnβ ( x)),
                                                      2
we have x, r2 and ru multiples of u2 whereas u4 is not. Otherwise, Lemma 17
makes it possible to conclude with a similar reasoning.
                                           Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
140      ·      Pierre Roux

5.    DIVISION
5.1    Even Radix
Figueroa [Fig95] proved double rounding to be innocuous for the division if it is
performed with a precision at least twice as large. Although his proof only addressed
radix 2, this holds for any even radix.
   Theorem 29. [Fig95] (double round div FLX) When β is even, for p1 , p2 ∈ Z,
if p2 ≥ 2p1 , then for ◦1 and ◦2 roundings to nearest, with any tie, respectively in
FLXp1 and FLXp2 :
                      ∀x, y ∈ FLXp1 , y 6= 0 ⇒ ◦1(◦2(x/y)) = ◦1(x/y) .
     This applies to IEEE754 binary32/64 when no underflow occurs.
   Remark 30. [Fig95] This bound 2p1 is optimal as shown by the following coun-
terexample. With radix β = 2, for p1 = 4, p2 = 7, we have for x = 1.0002 and
y = 0.11112
             ◦1(◦2(x/y)) = ◦1(◦2(1.000100010 . . .2 )) = ◦1(1.00010002 ) = 1.0002
which differs from ◦1(1.000100010 . . .2 ) = 1.0012 when ◦1 and ◦2 are roundings to
nearest with tie-break to even.
                                               ϕ2       ulp (x/y)
   Proof. Denoting midpϕ1(x/y) with m and       2     with u, from Lemma 16,
the result holds when x/y < m − u and similarly when x/y > m + u. When
m − u ≤ x/y < m, an argument similar to the one used in proof of Theorem 25
allows to conclude and similarly when m < x/y ≤ m + u. Finally, when x/y = m,
since β is even, x/y ∈ Fϕ2 as soon as ϕ2 (x/y) < ϕ1 (x/y) which is implied by
p2 ≥ 2p1 and p1 ≥ 1.
  The previous result still applies to IEEE754 binary32/64 even in case of underflow
(double round div FLT).
5.2    Odd Radix
The result does not hold for odd radices.
   Remark 31. The result does not hold for odd radices, even for arbitrarily large
p2 , as shown by the following counterexample. With radix β = 3, for p1 = 4, p2 ≥ 4,
p2 odd, for x = 20013 and y = 23
              ◦1(◦2(x/y)) = ◦1(◦2(1000.111 . . .3 )) = ◦1(1000.1 . . . 13 ) = 10003
which differs from ◦1(1000.111 . . .3 ) = 10013 when ◦1 and ◦2 are roundings to near-
est with tie-break to even.
  However, the result can be recovered by restricting the choice of tie-break rules
to directed rules such as round to nearest with tie-break away.
  Theorem 32. (double round div rna FLX) When β is odd, for p1 , p2 ∈ Z, if
p2 ≥ p1 , then for ◦A       A
                    1 and ◦2 roundings to nearest with tie-break away, respectively
in FLXp1 and FLXp2 :
                  ∀x, y ∈ FLXp1 , y 6= 0 ⇒ ◦A  A          A
                                                      
                                            1 ◦2 (x/y) = ◦1 (x/y) .

Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
                    Innocuous Double Rounding of Basic Arithmetic Operations               ·     141

Table I. Hypotheses under which double rounding is innocuous for the FLX format (β is the radix
and p1 and p2 are the two precisions involved (c.f., Section 1 for formal definitions)). Division for
odd radices is only valid with tie-break away.
                                    β=2            even β ≥ 4        odd β
                        +, −    p2 ≥ 2p1 + 1        p2 ≥ 2p1        p2   ≥ p1
                         ×        p2 ≥ 2p1          p2 ≥ 2p1        p2   ≥ p1
                         /        p2 ≥ 2p1          p2 ≥ 2p1        p2   ≥ p1
                         √
                                p2 ≥ 2p1 + 2      p2 ≥ 2p1 + 1      p2   ≥ p1

Table II. Hypotheses under which double rounding is innocuous for the FLT format, in addition
to hypotheses of Table I (emin1 and emin2 are the two minimal exponents involved). Division for
odd radices is only valid with tie-break away.
                         β=2                            even β ≥ 4                   odd β
     +, −          emin2 ≤ emin1                      emin2 ≤ emin1              emin2   ≤ emin1
      ×           emin2 ≤ 2emin1                     emin2 ≤ 2emin1              emin2   ≤ emin1
      /        emin2 ≤ emin1 − p1 − 2            emin2 ≤ emin1 − p1 − 2          emin2   ≤ emin1
      √
            (emin2 ≤ emin1 − p1 − 2             (emin2 < emin1 − p1              emin2   ≤ emin1
             ∨ 2emin2 ≤ emin1 − 4p1 − 2)         ∨ 2emin2 ≤ emin1 − 4p1 )

Table III. Hypotheses under which double rounding is innocuous for the FTZ format, in addition
to hypotheses of Table I. Division for odd radices is only valid with tie-break away.
                                   even β                                odd β
            +, −         emin2 + p2 ≤ emin1 + 1             emin2 + p2    ≤ emin1 + p1
             ×          emin2 + p2 ≤ 2emin1 + p1            emin2 + p2    ≤ emin1 + p1
             /             emin2 + p2 < emin1               emin2 + p2    ≤ emin1 + p1
             √
                    2 (emin2 + p2 ) ≤ emin1 + p1 ≤ 1        emin2 + p2    ≤ emin1 + p1

   Proof. According to Lemma 18, the result holds when x/y is not a midpoint.
It also holds when it is, since the directed tie-break rule ensures that midpoints are
always rounded in the same direction.

6.    CONCLUSION
Although all the results were illustrated throughout the paper on the FLX floating-
point format (no underflows), they remain valid under mild conditions on the FLT
(gradual underflow) and FTZ (abrupt underflow) formats with underflows. Indeed,
all our Coq proofs are carried out at the level of Fϕ generic formats and results
for the FLX, FLT and FTZ formats with or without underflows are just corollaries
immediatly derived from the former. The precise hypotheses under which double
rounding has been proved innocuous are summarized in Tables I, II and III.
   Once the proofs, originally expressed in terms of bit patterns [Fig95], have been
translated to the more algebraic form seen throughout the paper, Coq proofs ap-
peared rather easy to perform. A key factor in this easyness was the high genericity
of the definitions in the Flocq library [BM11]. For instance, most of our proofs are
done for any rounding to nearest, for any tie, by just using the definition of round-
ing to nearest available in Flocq. It should also be noticed that the library appeared
pretty complete since no definition and only a few (13) lemmas were added to it.
                                                 Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
142     ·      Pierre Roux

   All our Coq developments (6.9 kloc) are part of the Flocq library [BM11] which
is available1 under an open source license at http://flocq.gforge.inria.fr/.
   Thus, we have not only formally proved previously-known results on innocuous-
ness of double rounding of basic arithmetic operations [Fig95] but also extended
them to radices other than two (such as radix 10 introduced in the last revision of
the IEEE754 norm [IEE08]) and underflow cases (either gradual or abrupt). To the
extent of author knowledge, although they belong to the folklore in the floating-
point arithmetic community, no proof, neither pen and paper nor mechanical, can be
found in the literature. In particular, it is now formally proved that double round-
ing of addition/subtraction, multiplication, division, and square root, is innocuous
for IEEE754 binary32 format when using binary64 as intermediate precision. This
could, for instance, be used in the verified C compiler Compcert [BJLM13] to per-
form constant folding while cross compiling for a binary32 target on a binary64
host. This may also be of interest for some processors, such as the one using the
POWER (ancestor of PowerPC) instruction set, which have binary64 arithmetic
instructions but lack binary32 ones. The RAD6000, a radiation-hardened board
used onboard various spacecrafts, belongs to this category.

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

References
[BCHPM04] Yves Bertot, Pierre Castéran, Gérard (informaticien) Huet, and Chris-
          tine Paulin-Mohring. Interactive theorem proving and program devel-
          opment : Coq’Art : the calculus of inductive constructions. Texts
          in theoretical computer science. Springer, Berlin, New York, 2004.
          Données complémentaires http://coq.inria.fr.
[BJLM13]  Sylvie Boldo, Jacques-Henri Jourdan, Xavier Leroy, and Guillaume
          Melquiond. A Formally-Verified C Compiler Supporting Floating-
          Point Arithmetic. In Alberto Nannarelli, Peter-Michael Seidel, and
          Ping Tak Peter Tang, editors, IEEE Symposium on Computer Arith-
          metic, pages 107–115. IEEE Computer Society, 2013.
[BM11]    Sylvie Boldo and Guillaume Melquiond. Flocq: A Unified Library for
          Proving Floating-point Algorithms in Coq. In Proceedings of the 20th
          IEEE Symposium on Computer Arithmetic, pages 243–252, Tübingen,
          Germany, July 2011.
[Coq12]   The Coq development team. The Coq proof assistant reference man-
          ual, 2012. Version 8.4.
[Fig95]   Samuel A. Figueroa. When is Double Rounding Innocuous? SIGNUM
          Newsl., 30(3):21–26, July 1995.
[IEE08]   IEEE Computer Society. IEEE Standard for Floating-Point Arith-
          metic. IEEE Standard 754-2008, 2008.
[MDMM13] Érik Martin-Dorel, Guillaume Melquiond, and Jean-Michel Muller.
          Some issues related to double rounding. BIT Numerical Mathematics,
          53(4):897–924, 2013.
Journal of Formalized Reasoning Vol. 7, No. 1, 2014.
