# Computing the exact sign of sums of products with floating point arithmetic

Walter F. Mascarenhas  Departamento de Computação, IME
Universidade de São Paulo, Brazil

###### Abstract

In computational geometry, the construction of essential primitives like
convex hulls, Voronoi diagrams and Delaunay triangulations require the
evaluation of the signs of determinants, which are sums of products. The same
signs are needed for the exact solution of linear programming problems and
systems of linear inequalities. Computing these signs exactly with inexact
floating point arithmetic is challenging, and we present yet another algorithm
for this task. Our algorithm is efficient and uses only of floating point
arithmetic, which is much faster than exact arithmetic. We prove that the
algorithm is correct and provide efficient and tested C++ code for it.

##  1 Introduction

From a computational geometry perspective, we address the following problem:
given a point $C$, defined explicitly or implicitly as the intersection of two
lines, we want to decide on which side of a line $AB$ it lies using only
floating point arithmetic (or to find whether it is on the line.) This problem
is classic. It plays a fundamental role in the computation of convex hulls,
Voronoi diagrams and Delaunay triangulations [1, 14, 16]. The essence of the
geometric issues is described in Figure 1.

$C$$A$$B$rounded $C$ Figure 1: Rounding errors may lead to the conclusion that
the point $C$ lies in the wrong side of the line $AB$.

From an algebraic perspective, the determination of the position of the point
with respect to the line is equivalent to the computation of the sign of a sum
of products

$S=\sum_{i=1}^{n}\prod_{j=1}^{n_{i}}a_{ij}.$

The inexactness of floating point arithmetic makes it hard to compute this
sign exactly in some cases. There are already several references about the
computation sums of floating point numbers [3, 4, 5, 6, 7, 8, 10, 12, 13, 15,
16, 17], with applications in computations geometry. The idea of decomposing
products in sums is also as old as [2, 11], and the fourth chapter of [14]
presents an algorithm for obtaining the signs of such sums. However, all
theses references differ in one way or another from the present article. In
summary, none of them gives a complete solution to the problem we solve here,
using IEEE floating point arithmetic and taking underflow into account. We
present complete algorithm for this task, in gory detail. In particular, we
take as input the factors $a_{ij}$ of the products $p_{i}$ and before applying
most of the other algorithms we would need to obtain an exact representation
of the products $p_{i}$, or assume that such a product does not underflow. The
algorithm was implemented in C++, and it was carefully tested, and is
supported by a detailed theory. It is quite likely that the theoretical
results presented here can be generalized to bases other than two and other
rounding modes. We are not interested in such generalizations. The goal of our
theory is only to provide a justification to what matter most us: the
practical algorithms we present, and our theoretical results are only meant to
justify our algorithms.

Finite floating point numbers suffice to solve our problem in the situations
usually found in practice, and we only care about such numbers, which are
elements of a set ${\cal{F}}$. Deviating from the tradition, we assume that
the arithmetic operations $\mathrm{op}\in\left\\{{+,-,*}\right\\}$ are
performed rounding up. For instance, we assume that the subtraction of
floating point numbers is defined as

$b\ominus a={\mathrm{up}}\\!\left(b-a\right)$

where $\mathrm{up}$ is the function from
$\left({\min{\cal{F}},max{\cal{F}}}\right)\rightarrow{\cal{F}}$ defined by

${\mathrm{up}}\\!\left(x\right):=\min\left\\{{y\in{\cal{F}}\ \mathrm{with}\
y\geq x}\right\\},$ (1)

and the operations $\oplus$ and $\otimes$ are similar (we have no need for
expensive $\oslash$’s.) For brevity, we leave the result of the arithmetic
operation $a\ \mathrm{op}\ b$ undefined when

$a\ \mathrm{op}\ b\not\in\left({\min{\cal{F}},max{\cal{F}}}\right).$

In practice, when the rounding mode is not upwards already, enforcing our
assumption requires a function call to change the rounding mode upwards in the
beginning of the use our functions, and another function call to restore it
when we are done, to be polite. The cost of this two calls is amortized, but
cases like this make us believe that always rounding upwards (or always
rounding downards) is a good option for code requiring exact results. Our
Moore library [9] is a good example of this. Since it works with rounding mode
upwards by default, there would be no need for changes in the rounding mode
when executing the algorithms described here.

The motivation for this article is to have a bullet proof algorithm for
computing the sign of $S=\sum_{i=1}^{n}\prod_{j=1}^{n_{i}}a_{ij}$ with
floating point arithmetic. Such an algorithm is needed because sometimes
rounding errors may lead to wrong conclusions by naive algorithms. However,
floating point arithmetic is an excellent tool: such naive algorithms will be
correct most of the time, and we should care about both precision and
performance. For this reason, we believe that we should proceed in two steps:
we first try to compute the sign using a quick algorithm, which is unable to
compute the sign only in rare cases. In such rare cases we resort to a robust
algorithm, which can handle all cases but is more expensive. We present both
algorithms below. In Section 2 we describe a simple algorithm, which is quite
efficient but may not find the sign in rare cases. In Section 3 we present a
robust algorithm, which is more expensive and should be used only after a
quick algorithm was unable to find the sign.

Regarding the efficiency of our algorithms, we emphasize that the quick one
will suffice in the overwhelming majority of the cases found in practice, and
the robust algorithm will be an extreme safety measure. As a result, usually
the cost of evaluating the sign will be twice the cost of evaluating the sum
using naive floating point arithmetic, plus the cost of
$\sum_{i=1}^{n}n_{i}-n$ branches, plus the change of rounding modes when they
are necessary. This is a ${O}\\!\left(\sum_{i=1}^{n}n_{i}\right)$ cost, with a
small constant. However, in rare cases in which the robust algorithm is
necessary, the cost can grow exponentially with the $n_{i}$.

Finally, the actual code is implemented using the template features of C++,
and exploiting the details of this language it is possible to generate code
that will only resort to more expensive operations in the rare situations in
which they are needed.

##  2 Quick sign

In this section we present a fast algorithm which finds the sign of the sum of
products $S:=\sum_{i=1}^{n}\prod_{j=1}^{n_{i}}a_{ij}$ in most cases, but which
may be inconclusive sometimes. The algorithm returns a sign
$s\in\left\\{{-1,0,1}\right\\}$. If $s\neq 0$ then it definitely is the sign
of the sum. However, if $s=0$ then the sign can be anything. In this case we
must resort to a more expensive algorithm to find the sign, like the one in
the next section.

As in the rest of this article, we use floating point arithmetic rounding
upwards, with at most two changes in rounding mode. In essence, for each
product $p_{i}=\prod_{j=1}^{n_{i}}a_{ij}$ the algorithm computes numbers
$d_{i}$ and $u_{i}$ such that

$-d_{i}\leq p\leq u_{i}.$

If $\sum d_{i}<0$ then $S$ is positive and if $\sum u_{i}<0$ then $S$ is
negative, otherwise we cannot decide, and the algorithm returns $0$. The
algorithm is described as Algorithm 2 below, which uses the auxiliary
Algorithm 1.

Algorithm 1 Auxiliary Algorithm for bounding aij * x…

procedure quick_prod(d, u, aij, x…)

if aij $\leq$ 0 then $\triangleright$ $-d\leq\prod_{k=0}^{j-1}a_{ik}\leq u$

if aij $=$ 0 then

return 0,0

else$\triangleright$ $-u\leq-\prod_{k=0}^{j-1}a_{ik}\leq d$

return quick_prod(u, d, -aij, x…) $\triangleright$
$a_{ij}\prod_{k=0}^{j-1}a_{ik}=-a_{ij}*-\prod_{k=0}^{j-1}a_{ik}$

d $\leftarrow$ d * aij $\triangleright$ $-\left({d\otimes
a_{ij}}\right)\leq-\left({da_{ij}}\right)\leq\prod_{k=0}^{j}a_{ik}$

u $\leftarrow$ u * aij $\triangleright$ $\prod_{k=0}^{j}a_{ik}\leq ua_{ij}\leq
u\otimes a_{ij}$

if is_empty(x) then

return d, u

else

return quick_prod(d, u, x…)

end procedure

Algorithm 2 A quick algorithm for computing the sign of Sum p…

procedure quick_sign(p…)

round(up) $\triangleright$ round object, resets the rounding mode on exit

su $\leftarrow$ 0 $\triangleright$ su is an upper bound on $S$

sd $\leftarrow$ 0 $\triangleright$ sd is an upper bound on $-S$.

for  each p do

d, u $\leftarrow$ quick_prod(-1, 1, p)

su $\leftarrow$ su + u

sd $\leftarrow$ sd + d

return sd $<0$ ? 1 : (su $<$ 0 ? -1: 0)

end procedure

##  3 Robust sign

This section describes Algorithm 3, which computes the sign of
$S=\sum_{i=1}^{n}\prod_{j=1}^{n_{i}}a_{ij}$, using binary floating point
arithmetics which have subnormal numbers. The most relevant arithmetics in
this class are the ones covered by he IEEE 754 and IEEE 854 standards. The
algorithm is presented in the last page of the article. It assumes that there
is not overflow in the products $p_{i}=\sum_{j=1}^{n_{i}}a_{ij}$, but
underflow is handled correctly. It also assumes that

$\sum_{i\ \mathrm{with}\ n_{i}=1}1+\sum_{i\ \mathrm{with}\
n_{i}>1}2^{n_{i}-2}<1/\epsilon,$ (2)

where $\epsilon$ is the machine precision. In practice, the largest $\epsilon$
we care about corresponds to the type float. In this case $1/\epsilon=2^{23}$
and the algorithm could be used to compute the signs of determinants of
dimensions up to 8, because if $n=8!$ and $n_{i}=8$ then

$\sum_{i=1}^{n}2^{n_{i}-2}=8!\times 2^{6}=2580480<8388608=2^{23}.$

Since $n_{i}=8$ is already more than enough for the usual applications in
computational geometry, we have no motive to make the algorithm more
complicated than it already is in order to relax the condition
$\left({\ref{boundni}}\right)$.

The algorithm is based upon two lemmas. Lemma 1 is about the exact computation
of the difference $b-a$ of floating point numbers. There are versions of this
lemma since the late 1960s [2, 7, 8], but we prove it here for the case in
which we round upwards for completeness, and because the details are not
obvious (as stated, Lemma 1 is false if we round downwards for instance.) In
essence, it states that we can represent $b-a$ exactly as the difference of
two floating point numbers $c$ and $e$, with the additional feature that $e$
is much smaller than $c$. As a result, in most cases we can base our decision
regarding signs on $c$, and $e$ is used only in the rare cases in which
knowing $c$ is not sufficient.

###### Lemma 1

if $a,b\in{\cal{F}}$, $0<a<b$ and

$c:=b\ominus a,\hskip 28.45274ptd:=b\ominus c\hskip
28.45274pt\mathrm{and}\hskip 28.45274pte:=a\ominus d$

then

$b-a=c-e\hskip 28.45274pt\mathrm{and}\hskip 28.45274pt0\leq e<c\epsilon,$ (3)

where $\epsilon$ is the machine precision. $\blacktriangle$

Lemma 2 is the analogous to Lemma 1 for multiplication, but it is more subtle.
It relies on the fused multiply add operation ($\mathrm{fma}$), which is
available in most processors and programming language these days. In other
words, we assume that given $a,x,y\in{\cal{F}}$ such that
$ax+y\in\left({\min{\cal{F}},\max{\cal{F}}}\right)$ we can compute

${\mathrm{fma}}\\!\left(a,x,y\right):={\mathrm{up}}\\!\left(ax+y\right).$ (4)

It is well known that using the fma we can represent the product $ab$ as the
difference of two floating point numbers, but we are not aware of proofs (or
even statements) of results describing conditions under which this
representation is exact when rounding upwards. It is important to notice that
such conditions are necessary, because Lemma 2 may not hold if the condition
$\left({\ref{hypo_fma}}\right)$ is violated. In order to state the
decomposition result for multiplications we need to define some constants that
characterize the floating point arithmetic:

  * •

$\nu$ is the smallest positive normal element of ${\cal{F}}$ and $\nu\epsilon$
is the smallest positive element of ${\cal{F}}$

  * •

$\epsilon\in{\cal{F}}$ is the machine precision, that is, $1+\epsilon$ is the
successor of $1$ in ${\cal{F}}$.

  * •

$\sigma$ is the largest power of two in ${\cal{F}}$, and we assume that
$1/\sigma\in{\cal{F}}$.

Using the constants above we define the threshold

$\tau:=2\nu/\epsilon.$ (5)

The values of these constants for the usual arithmetics are presented in Table
1. By inspecting this table, readers will notice that the following
assumptions used in Lemma 2 are satisfied:

$\sigma\epsilon^{2}\geq 2\hskip 28.45274pt\mathrm{and}\hskip
28.45274pt\sigma^{2}\nu\epsilon>2.$ (6)

| float | double | long double | quad
---|---|---|---|---
$\nu$ | 1.2e-38 | 2.2e-308 | 3.4e-4932 | 3.3e-4932
$\epsilon$ | 1.2e-07 | 2.2e-016 | 1.1e-0019 | 1.9e-0034
$\sigma$ | 1.7e+38 | 9.0e+307 | 5.9e+4931 | 5.9e+4931
$\tau$ | 2.0e-31 | 2.0e-292 | 6.2e-4913 | 3.4e-4898
Table 1: Values of the constants for the usual arithmetics.

We now state the decomposition lemma for multiplications.

###### Lemma 2

Consider $a,b\in{\cal{F}}$, with $0<a\leq b$, for which $a\otimes b$ is
defined. If

$c:=a\otimes b>\tau,$ (7)

then

$ab=c-d\hskip 28.45274pt\mathrm{for}\hskip
28.45274ptd:={\mathrm{fma}}\\!\left(-a,b,c\right)\geq 0.$ (8)

If the condition $\left({\ref{hypo_fma}}\right)$ is not satisfied then
$\tilde{a}:=\sigma\otimes a=\sigma a\in{\cal{F}}$ and if

$\tilde{c}:=\tilde{a}\otimes b>\tau,$ (9)

then

$ab=\sigma^{-1}\tilde{c}-\sigma^{-1}\tilde{d}\hskip
28.45274pt\mathrm{for}\hskip
28.45274pt\tilde{d}:={\mathrm{fma}}\\!\left(-\tilde{a},b,\tilde{c}\right)\geq
0.$ (10)

Finally, if Equations $\left({\ref{hypo_fma}}\right)$ and
$\left({\ref{hypo_fma2}}\right)$ are not satisfied and the first inequality in
Equation $\left({\ref{asa}}\right)$ holds then $\hat{b}:=\sigma\otimes
b=\sigma b\in{\cal{F}}$, and if the second inequality in Equation
$\left({\ref{asa}}\right)$ holds then

$ab=\sigma^{-2}\hat{c}-\sigma^{-2}\hat{d}\hskip 14.22636pt\mathrm{for}\hskip
14.22636pt\hat{c}:=\tilde{a}\otimes\hat{b}\hskip 14.22636pt\mathrm{and}\hskip
14.22636pt\hat{d}:={\mathrm{fma}}\\!\left(-\tilde{a},\hat{b},\hat{c}\right)\geq
0.$ (11)

In summary, if $0<a\leq b$ and $a\otimes b$ is defined then there exists
$e\in\left\\{{0,1,2}\right\\}$ and $c,d\in{\cal{F}}$ such that

$ab=\sigma^{-e}c-\sigma^{-e}d.$ (12)

$\blacktriangle$

In words, Lemma 2 shows that we may fail twice in trying to represent exactly
the product $ab$ as a difference of two floating point numbers, but the third
time is a charm: we finally can represent $ab$ exactly as the scaled
difference of two floating point numbers. Scaling is essential here in order
to deal with underflow. We scale numbers by multiplying them by the constant
$\sigma$, which is a power of two. Such scaling does not introduce rounding
errors, but requires some book keeping. In C++ we can keep the books using an
struct like

    struct scaled_number {
      scaled_number(T t, int exp)
      T t;
      int exp;
    };

where T is the type of the floating point numbers. An scaled number s
represents

$x={\mathrm{value}}\\!\left(\texttt{s}\right)=\sigma^{-\texttt{s.exp}}\
\texttt{s.t}.$

We keep the scaled numbers in two heaps, one for the positive values, called
${\cal{P}}$, and another for the negative values, called ${\cal{N}}$. In the
${\cal{N}}$ heap we store the absolute value of the corresponding numbers, so
that the t field in our scaled numbers is always positive, and
$\texttt{exp}\leq 0$. The elements in the heaps are sorted in increasing order
according to the following comparison function:

    bool is_less(scaled_number x, scaled_number y) {
      if( x.exp > y.exp ) return true;
      if( x.exp < y.exp ) return false;Ψ
      return x.t < y.t;
    }

In order to ensure the consistency of the order above we only push two kinds
of scaled numbers in our heaps:

$\displaystyle\mathrm{Scaled\ numbers}\ \texttt{s}\ \mathrm{with}$
$\displaystyle\texttt{s.exp}=0\ \mathrm{and}\ \tau<\texttt{s.t},$ (13)
$\displaystyle\mathrm{Scaled\ numbers}\ \texttt{s}\ \mathrm{with}$
$\displaystyle\ \texttt{s.exp}>0\ \mathrm{and}\
\tau<\texttt{s.t}\leq\sigma\tau.$ (14)

We assume that these conditions are enforced by the constructor of
scaled_numbers, which is only called with positive $t$’s. We then have the
following Lemma

###### Lemma 3

Under the conditions $\left({\ref{scale0}}\right)$ and
$\left({\ref{scale1}}\right)$ for scaled numbers x and y we have

${\mathrm{value}}\\!\left(\texttt{x}\right)<{\mathrm{value}}\\!\left(\texttt{y}\right)\
\Leftrightarrow\ {\texttt{is\\_less}}\\!\left(\texttt{x},\texttt{y}\right).$

$\blacktriangle$

We have now all the ingredients to describe our algorithm. It uses an
auxiliary function split which splits each product
$p_{i}=\prod_{j=1}^{n_{i}}a_{ij}$ as a sum of $2^{n_{i}-1}$ scaled numbers
using Lemma 2 (given this lemma, writing such a function is trivial.) If
$n_{i}>1$ then half of the parts in which $p_{i}$ is split will be negative
and the other half will be positive. Therefore, if $n_{i}>1$ then $p_{i}$
contributes $2^{n_{i}-2}$ scaled numbers to each heap. As a result, the left
hand side of Equation $\left({\ref{boundni}}\right)$ is the maximum number of
elements which we will have on each heap, and the condition
$\left({\ref{boundni}}\right)$ ensures that this number does not exceed
$1/\epsilon$.

Once the heaps are filled with products, we start to compare them. While

$s_{n}:={\mathrm{size}}\\!\left({\cal{N}}\right)>0\hskip
28.45274pt\mathrm{and}\hskip
28.45274pts_{p}:={\mathrm{size}}\\!\left({\cal{P}}\right)>0$

we pop the top elements $n$ and $p$ of ${\cal{N}}$ and ${\cal{P}}$ and compare
them. If $p>s_{n}*n$ then the sign of the sum is $1$. If $n>s_{p}*p$ then the
sign of the sum is $-1$, otherwise, conceptually, we use Lemma 1 to split
$p-n$ and push the parts back into the heaps. There is a catch in this
argument in that $n$ and $p$ are scaled numbers, which may have different
exponents. If these exponents differ by more than one then the numbers on the
heap with the numbers with the largest exponent are negligible and we are
done. When the exponents differ by one we multiply the t field of the one with
the largest exponent by $\sigma$, reducing both numbers to the same exponent.
This multiplication by $1/\sigma$ may be inexact, but this inexactness is
harmless. For instance, when $p$ has the largest exponent and the
multiplication $z:=p.t\leftarrow p.t\otimes 1/\sigma$ is inexact, then Lemma 4
in Section 4 implies that $z\leq\nu$ and Equations
$\left({\ref{scale0}}\right)$ and $\left({\ref{scale1}}\right)$ yield

$n.t>\tau=\frac{2}{\epsilon}\nu\geq\frac{2}{\epsilon}z\geq 2s_{p}z,$ (15)

and we will reach the correct conclusion that the sign of the sum is $-1$ even
if we use the incorrect $z$. Once we have both numbers with the same exponent,
simply split the difference $p.t-n.t$ and adjust the exponents of the results
consistently.

Finally, the algorithm terminates because we have two possibilities after we
reduce $p$ and $n$ to the same exponent. When $n.t\leq p.t$ (and the case
$n.t>p.t$ is analogous.):

  * (i)

If $p.t\ominus p.n<p.t$ then the largest $t$ field decreases, and this can
only happen a finite number of times.

  * (ii)

If $p.t\ominus p.n=p.t$ then Lemma 1 implies that $e=p.n<\epsilon p.t$. Since
the number of elements in ${\cal{N}}$ is at most $1/\epsilon$ by bound
$\left({\ref{boundni}}\right)$, this implies that $p.t>s_{n}*p.n$, and the
algorithm returns $1$ due to this condition.

This is only a high level description of the algorithm. A reasonably detailed
version of it is presented in the last page of this article. The actual code
is a bit more involved, due to optimizations which replace scaled_numbers by
plain floating point numbers when possible. Readers interested in the
implementation details should look at the C++ code available as supplementary
material to the arxiv version of this article. This code is distributed under
the Mozilla Public License 2.0.

Algorithm 3 Sign of Sums of Products

1:procedure sign(p…) $\triangleright$ sign of the sum p[0] + p[1] +…

2: round r(up) $\triangleright$ Rounding object, resets the mode on exit

3: heap neg, pos

4:

5: for  each p do

6: split(neg, pos, p) $\triangleright$ Each p is a product

7:

8: while true do

9: sn, sp $\leftarrow$ size(neg), size(pos)

10:

11: if sn $=$ 0 then

12: return (sp $=$ 0) ? 0 : 1

13: if sp $=$ 0 then

14: return -1

15:

16: n, p $\leftarrow$ pop(neg), pop(pos)

17:

18: if n.exp $\leq$ p.exp then $\triangleright$ small n.exp $\Rightarrow$
large n $=2^{-\mathrm{n.exp}}$ n.t

19: if n.exp $<$ p.exp then

20: if n.exp $<$ p.exp - 1 then

21: return -1 $\triangleright$ p is too small

22: p.t $\leftarrow$ p.t * $1/\sigma$ $\triangleright$ $p.t/\sigma$ is exact
or -1 is returned in line 24

23: if n.t $>$ sp * p.t then

24: return -1

25:

26: if n.t $>$ p.t then

27: c, e $\leftarrow$ lemma_1(p.t, n.t) $\triangleright$ split n.t - p.t as in
Lemma 1

28: neg $\leftarrow$ scaled_number(c, n.exp)

29: if e $>$ 0 then

30: pos $\leftarrow$ scaled_number(e, n.exp)

31: else

32: if p.t $>$ n.t then

33: c, e $\leftarrow$ lemma_1(n.t, p.t) $\triangleright$ split p.t - n.t as in
Lemma 1

34: pos $\leftarrow$ scaled_number(c, n.exp)

35: if e $>$ 0 then

36: neg $\leftarrow$ scaled_number(e, n.exp)

37: else

38: if n.exp $\geq$ p.exp $+$ 1 then

39: return 1 $\triangleright$ n is too small

40: n.t $\leftarrow$ n.t * 1/$\sigma$ $\triangleright$ $n.t/\sigma$ is exact
or 1 is returned in line 42

41: if p.t $>$ sn * n.t then

42: return 1

43: c, e $\leftarrow$ lemma_1(n.t, p.t) $\triangleright$ split p.t - n.t as in
Lemma 1

44: pos $\leftarrow$ scaled_number(c, p.exp)

45: if e $>$ 0 then

46: neg $\leftarrow$ scaled_number(e, p.exp)

47: end while

48:end procedure

##  4 Proofs

Here we prove the results stated above and two auxiliary results. Our proofs
use the following characteristics shared by the usual binary floating point
arithmetics with subnormal numbers, like the ones covered by the IEEE
standards 754 and 854. There are three kinds of elements in the set
${\cal{F}}$ of finite floating point numbers:

  * •

$0$ is a floating point number.

  * •

$x\in{\cal{F}}$ if and only if $-x\in{\cal{F}}$.

  * •

The subnormal numbers $x$ have absolute value of the form

$x=\nu\epsilon m,\hskip 28.45274pt\mathrm{with}\hskip
28.45274ptm\in{\mathds{N}}{}\hskip 28.45274pt\mathrm{and}\hskip
28.45274pt1\leq m<1/\epsilon.$ (16)

  * •

The normal numbers $x$ have absolute value of the form

$x=2^{e}\nu\epsilon m$ (17)

for integers $e$ and $m$ such that

$0\leq e\leq
e_{\max}:={\log_{2}}\\!\left(\sigma\right)-{\log_{2}}\\!\left(\nu\right)\hskip
28.45274pt\mathrm{and}\hskip 28.45274pt1\leq\epsilon m<2.$ (18)

We use two auxiliary lemmas, and the proofs are presented after the statements
of these lemmas. The lemmas are proved in the order in which they were stated.

###### Lemma 4

If $x\in{\cal{F}}$ and $2{\left|{x}\right|}\leq\max{\cal{F}}$ then
$2x\in{\cal{F}}$. Therefore, if $e$ is a positive integer and
$2^{e}{\left|{x}\right|}\leq\max{\cal{F}}$ then $2^{e}x\in{\cal{F}}$.
Similary, if $k\geq 0$ and $2^{-k}{\left|{x}\right|}\geq\nu$ then
$2^{-k}\otimes x=2^{-k}x$ is normal. $\blacktriangle$

###### Lemma 5

If the integer $\ell\geq 0$ and the real number $x$ are such that

$2^{\ell}\leq\frac{x}{\nu}\leq 2^{\ell+1}\leq\sigma$ (19)

then

${\mathrm{up}}\\!\left(x\right)=2^{\ell}\nu\epsilon\left\lceil{{\frac{x}{2^{\ell}\nu\epsilon}}}\right\rceil.$
(20)

$\blacktriangle$

Proof of Lemma 1. If $b=\nu\epsilon m_{b}$ is either subnormal or normal with
a minimum exponent, then $a$ is also of the form $a=\nu\epsilon m_{a}$,

$c=\nu\epsilon\left({m_{b}-m_{a}}\right),\hskip 28.45274ptd=\nu\epsilon
m_{a},\hskip 28.45274pt\ e=0\hskip 28.45274pt\mathrm{and}\hskip
28.45274ptb-a=c=c+e,$

and Equation $\left({\ref{thesis_sub}}\right)$ holds. We can then assume that

$b=2^{e_{b}}\nu\epsilon m_{b}\hskip 28.45274pt\mathrm{with}\hskip
28.45274pte_{b}>0\hskip 28.45274pt\mathrm{and}\hskip 28.45274pt1\leq\epsilon
m_{b}<2,$ (21)

and

$a=2^{e_{a}}\nu\epsilon m_{a}\hskip 28.45274pt\mathrm{with}\hskip
28.45274pte_{a}\geq 0\hskip 28.45274pt\mathrm{and}\hskip 28.45274pt1\leq
m_{a}<2/\epsilon.$

It follows that

$b-a=\nu\epsilon h\hskip 28.45274pt\mathrm{for}\hskip
28.45274pth:=2^{e_{b}}m_{b}-2^{e_{a}}m_{a}>0.$ (22)

If $h<1/\epsilon$ then $b-a$ is subnormal and

$c=b\ominus a=b-a\ \Rightarrow\ d=a\ \Rightarrow\ e=0\ \Rightarrow\ b-a=c-e,$

Equation $\left({\ref{thesis_sub}}\right)$ holds and we are done. Let us then
assume that $h\geq\epsilon$ and let $\ell\geq 0$ be the integer such that

$2^{\ell}\leq\frac{b-a}{\nu}=\epsilon h\leq 2^{\ell+1}.$ (23)

By Lemma 5,

$c=b\ominus
a=2^{\ell}\nu\epsilon\left\lceil{{2^{-\ell}h}}\right\rceil=2^{\ell}\nu\epsilon\overline{c}$
(24)

for

$\overline{c}:=\left\lceil{{2^{e_{b}-\ell}m_{b}-2^{e_{a}-\ell}m_{a}}}\right\rceil.$
(25)

Equation $\left({\ref{squeeze_sub}}\right)$ leads to the bound

$h=2^{e_{b}}m_{b}-2^{e_{a}}m_{a}\geq 2^{\ell}/\epsilon$ (26)

and Equation $\left({\ref{hsub}}\right)$ shows that

$2^{e_{b}}m_{b}\geq 2^{\ell}/\epsilon\ \Rightarrow\
2^{e_{b}}\geq\frac{2^{\ell}}{\epsilon m_{b}}>2^{\ell-1}\Rightarrow
e_{b}\geq\ell,$

and $2^{e_{b}-\ell}m_{b}$ is integer.

If $2^{e_{a}-\ell}m_{a}<1$ then

$\overline{c}=2^{e_{b}-\ell}m_{b},\hskip 14.22636ptc=b,\hskip
14.22636ptd=b\ominus c=0,\hskip 14.22636pte=a\ominus d=a\ \ \Rightarrow\ \
b-a=c-e,$

and the first part of Equation $\left({\ref{thesis_sub}}\right)$ holds.
Moreover, Equation $\left({\ref{hsub}}\right)$ yields

$e=a=2^{e_{a}}m_{a}\nu\epsilon<2^{\ell}\nu\epsilon<\left({2^{e_{b}}m_{b}\epsilon}\right)\nu\epsilon=\epsilon
b=\epsilon c,$

and the second part of Equation $\left({\ref{thesis_sub}}\right)$ also holds.
Therefore, we can assume that

$2^{\ell-e_{a}}\leq m_{a}<2/\epsilon.$ (27)

If $e_{a}\geq\ell$ then $2^{e_{a}-\ell}$ is integer,

$\overline{c}=2^{e_{b}-\ell}m_{b}-2^{e_{a}-\ell}m_{a},\hskip
14.22636ptc=b-a,\hskip 14.22636ptd=b\ominus c=a,\hskip 14.22636pte=a\ominus
d=0\ \ \Rightarrow\ \ b-a=c-e$

and Equation $\left({\ref{thesis_sub}}\right)$ holds again. Therefore, we can
assume that $e_{a}<\ell$ and

$q:=\ell-e_{a}>0.$ (28)

Equation $\left({\ref{bma}}\right)$ shows that $2^{-q}m_{a}\geq
2^{e_{a}-\ell}\times 2^{\ell-e_{a}}=1$ and the integers

$\overline{a}:=\left\lfloor{{2^{-q}m_{a}}}\right\rfloor\hskip
28.45274pt\mathrm{and}\hskip 28.45274pt\underline{a}=m_{a}-2^{q}\overline{a}$

are such that

$m_{a}=2^{q}\overline{a}+\underline{a},\hskip 14.22636pt1\leq\overline{a}\leq
m_{a}/2<1/\epsilon\hskip 14.22636pt\mathrm{and}\hskip
14.22636pt0\leq\underline{a}<2^{q}\leq m_{a}<2/\epsilon.$ (29)

It follows that

$2^{e_{b}-\ell}m_{b}-2^{e_{a}-\ell}m_{a}=2^{e_{b}-\ell}m_{b}-2^{-q}m_{a}=2^{e_{b}-\ell}m_{b}-\overline{a}-2^{-q}\underline{a}$

and the bound $0\leq\underline{a}<2^{q}$ leads to

$\overline{c}=\left\lceil{{2^{e_{b}-\ell}m_{b}-2^{e_{a}-\ell}m_{a}}}\right\rceil=2^{e_{b}-\ell}m_{b}-\overline{a}$

and

$c=b\ominus a=b-2^{\ell}\nu\epsilon\overline{a}=b-\hat{a}\hskip
28.45274pt\mathrm{for}\hskip
28.45274pt\hat{a}:=2^{\ell}\nu\epsilon\overline{a}.$ (30)

The second bound in Equation $\left({\ref{abar}}\right)$ shows that
$\nu\epsilon\overline{a}$ is subnormal and Lemma 4 shows that
$\hat{a}\in{\cal{F}}$. This implies that

$d=b\ominus c=b-c=\hat{a}.$ (31)

Additionally,

$a-d=2^{e_{a}}\nu\epsilon
m_{a}-2^{\ell}\nu\epsilon\overline{a}=2^{e_{a}}\nu\epsilon\left({m_{a}-2^{q}\overline{a}}\right)=2^{e_{a}}\nu\epsilon\underline{a}=:\tilde{a}.$
(32)

If $\underline{a}=0$ then $\tilde{a}\in{\cal{F}}$. If
$\underline{a}<1/\epsilon$ then the same argument used for $\hat{a}$ shows
that $\tilde{a}\in{\cal{F}}$, and if $1/\epsilon\leq\underline{a}<2/\epsilon$
then $\nu\epsilon\underline{a}$ is normal and Lemma 4 shows that
$\tilde{a}\in{\cal{F}}$. Therefore, in all cases for $\underline{a}$ in
Equation $\left({\ref{abar}}\right)$ we have that $\tilde{a}\in{\cal{F}}$ and

$e=a\ominus d=a-d=\tilde{a}.$ (33)

Equations $\left({\ref{subd}}\right)$ and $\left({\ref{subad}}\right)$ show
that $\hat{a}+\tilde{a}=a$ and Equations $\left({\ref{subc}}\right)$ and
$\left({\ref{sube}}\right)$ yield

$c-e=b-\hat{a}-\tilde{a}=b-a,$

and the first part of Equation $\left({\ref{thesis_sub}}\right)$ holds.
Finally, Equations $\left({\ref{csub}}\right)$, $\left({\ref{cbar}}\right)$,
$\left({\ref{hsub}}\right)$,$\left({\ref{qsub}}\right)$ and
$\left({\ref{abar}}\right)$ imply that

$\epsilon c\geq
2^{\ell}\nu\epsilon=2^{q+e_{a}}\nu\epsilon>2^{e_{a}}\underline{a}\nu\epsilon=\tilde{a}.$

and the second part of Equation $\left({\ref{thesis_sub}}\right)$ holds.
$\blacksquare$

Proof of Lemma 2. Let us start with the case in which $c=a\otimes
b>\tau=2\nu/\epsilon$. In this case $b>\sqrt{2\nu/\epsilon}>\nu$ and $b$ is
normal. Therefore,

$b=2^{e_{b}}m_{b}\hskip 28.45274pt\mathrm{with}\hskip 28.45274pt1\leq\epsilon
m_{b}<2.$ (34)

$a$ can be normal or subnormal, but in both cases there exist integers
$e_{a}$, $\ell_{a}$ and $m_{a}$ with

$a=2^{e_{a}}m_{a},\hskip 28.45274pt2^{\ell_{a}}\leq m_{a}<2^{\ell_{a}+1}\hskip
28.45274pt\mathrm{and}\hskip 28.45274pt1\leq 2^{\ell_{a}}<1/\epsilon.$ (35)

If $\ell>0$ is the integer such that

$2^{\ell}\leq\frac{ab}{\nu}<2^{\ell+1}$ (36)

then

$2^{\ell+1}\geq\frac{ab}{\nu}>1/\epsilon\ \Rightarrow
2^{\ell}>1/\left({2\epsilon}\right)>1\ \ \Rightarrow\ell>0,$

and Lemma 5 and Equations $\left({\ref{eqb}}\right)$ and
$\left({\ref{eqa}}\right)$ show that

$a\otimes
b=2^{\ell}\nu\epsilon\left\lceil{{\frac{2^{e_{a}+e_{b}}m_{a}m_{b}}{2^{\ell}\nu\epsilon}}}\right\rceil=\left\lceil{{\frac{m_{a}m_{b}}{p}}}\right\rceil$
(37)

for

$p:=2^{\ell-e_{a}-e_{b}}\nu\epsilon.$ (38)

Since

$\frac{ab}{2^{\ell}v\epsilon}=\frac{m_{a}m_{b}}{p},$

Equation $\left({\ref{squeeze}}\right)$ implies that

$\frac{ab}{v}=\frac{2^{\ell}\epsilon m_{a}m_{b}}{p}<2^{\ell+1}$

and Equations $\left({\ref{eqb}}\right)$ and $\left({\ref{eqa}}\right)$ yield

$p>\epsilon m_{a}m_{b}/2=\left({m_{a}/2}\right)\left({\epsilon
m_{b}}\right)\geq 2^{\ell_{a}-1}\geq 1/2\Rightarrow p>1/2.$

Since $p$ is a power of two, this implies that $p\geq 1$ is integer. On the
other hand, Equation $\left({\ref{squeeze}}\right)$ implies that

$\frac{2^{\ell}\epsilon m_{a}m_{b}}{p}\geq 2^{\ell}$

and Equations $\left({\ref{eqb}}\right)$ and $\left({\ref{eqa}}\right)$ lead
to

$p\leq\epsilon m_{a}m_{b}<4/\epsilon\Rightarrow p\leq 2/\epsilon.$

Euclid’s division algorithm yields integers $q$ and $r$ such that

$m_{a}m_{b}=qp-r\hskip 28.45274pt\mathrm{and}\hskip 28.45274pt0\leq r<p\leq
2/\epsilon,$ (39)

and Equation $\left({\ref{sub6}}\right)$ shows that

$c=a\otimes
b=2^{\ell}\nu\epsilon\left\lceil{{q-p^{-1}r}}\right\rceil=2^{\ell}\nu\epsilon
q=ab+2^{\ell}p^{-1}\nu\epsilon r.$ (40)

Equations $\left({\ref{eqb}}\right)$ and $\left({\ref{eqa}}\right)$ imply that

$m_{a}m_{b}<2^{\ell_{a}+1}2/\epsilon<\frac{4}{\epsilon^{2}}\
\Rightarrow\frac{1}{m_{a}m_{b}}>\frac{\epsilon^{2}}{4},$

the bound $a\otimes b>2\nu/\epsilon$ yields $ab\epsilon>2\nu$ and

$2^{\ell}p^{-1}=\frac{2^{e_{a}+e_{b}}}{\nu\epsilon}=\frac{1}{\epsilon^{2}m_{a}m_{b}}\frac{2^{e_{a}+e_{b}}m_{a}m_{b}\epsilon}{\nu}=\frac{1}{\epsilon^{2}m_{a}m_{b}}\frac{ab\epsilon}{\nu}>\frac{2}{\epsilon^{2}m_{a}m_{b}}>\frac{1}{2}.$

Since $2^{\ell}p^{-1}$ is a power of 2, this implies that $2^{\ell}p^{-1}$ is
integer, and the last inequality in Equation $\left({\ref{mambf}}\right)$ and
Lemma 4 imply that

$d:=2^{\ell}p^{-1}\nu\epsilon r\in{\cal{F}}.$

Combining this with Equations $\left({\ref{otimes}}\right)$ we obtain that

$ab=2^{\ell}\nu\epsilon q-d=c-d$

and

$\left({-a}\right)b+c=d\in{\cal{F}}\Rightarrow{\mathrm{fma}}\\!\left(-a,b,c\right)=d.$

This completes the proof of Equation $\left({\ref{thesis_fma}}\right)$.

If $a\otimes b\leq\tau$ then

$a^{2}\leq ab\leq\tau\Rightarrow
a\leq\sqrt{\frac{2\nu}{\epsilon}}<1\Rightarrow\sigma
a<\sigma\Rightarrow\tilde{a}=\sigma a=\sigma\otimes a\in{\cal{F}},$

and using the same argument used to prove Equation
$\left({\ref{thesis_fma}}\right)$ we can prove Equation
$\left({\ref{thesis_fma2}}\right)$.

Finally, the smallest value possible for a positive floating point number is
$\nu\epsilon$, and if the conditions $\left({\ref{hypo_fma}}\right)$ and
$\left({\ref{hypo_fma2}}\right)$ are not satified then

$\tilde{a}=\sigma a\geq\sigma\nu\epsilon.$

The violation of condition $\left({\ref{hypo_fma2}}\right)$ and the first
inequality in Equation $\left({\ref{asa}}\right)$ yield

$b\leq\frac{2\nu}{\epsilon\tilde{a}}\leq\frac{2\nu}{\epsilon\sigma\nu\epsilon}\Rightarrow\sigma
b\leq\frac{2}{\epsilon^{2}}\leq\sigma.$

This bound implies that $\hat{b}=\sigma\otimes b=\sigma b\in{\cal{F}}$. We
also have that $\hat{b}\geq\sigma\nu\epsilon$ and the second inequality in
Equation $\left({\ref{asa}}\right)$ leads to

$\tilde{a}\hat{b}\geq\sigma^{2}\nu^{2}\epsilon^{2}=\frac{\sigma^{2}\nu\epsilon}{2}\frac{2\nu}{\epsilon}>\frac{2\nu}{\epsilon}=\tau.$

This condition allows us to use the same argument used to prove the validity
of Equations $\left({\ref{thesis_fma}}\right)$ and
$\left({\ref{thesis_fma2}}\right)$ to prove Equations
$\left({\ref{thesis_fma3}}\right)$, and this proof is complete. $\blacksquare$

Proof of Lemma 3. If $\texttt{x.exp}>\texttt{y.exp}$ then $\texttt{x.exp}>0$
and item (ii) implies that $\texttt{x.t}\leq\sigma\tau$. It follows that
$\texttt{x.exp}\geq\texttt{y.exp}+1$ and using item (i) we derive that

${\mathrm{value}}\\!\left(\texttt{x}\right)=\sigma^{-\texttt{x.exp}}\texttt{x.t}\leq\sigma^{-\texttt{y.exp}-1}\
\sigma\tau=\sigma^{-\texttt{y.exp}}\ \tau<\sigma^{-\texttt{y.exp}}\
\texttt{y.t}={\texttt{value}}\\!\left(y\right),$

and the function is_less returns the correct value in this branch. The branch
$\texttt{x.exp}<\texttt{y.exp}$ is analogous. Finally, in the last branch the
exponents cancel out and we are left with the correct comparison of the t
fields. $\blacksquare$

Proof of Lemma 4. If $x=0$ then $2x=0\in{\cal{F}}$. If $x$ is normal, then
Lemma 4 follows directly from Equation $\left({\ref{normal}}\right)$. If
$x=\nu\epsilon m$ is subnormal then there are two possibilities: If
${\left|{2m}\right|}<1/\epsilon$ then $2x=\nu\epsilon\left({2m}\right)$ is
also subnormal. If ${\left|{2m}\right|}\geq 1/\epsilon$ then
$2x=\nu\epsilon\left({2m}\right)$ is normal, because
$1/\epsilon\leq{\left|{m}\right|}<2/\epsilon$. The part of the Lemma regarding
division follows directly from the definition of normal number. $\blacksquare$

Proof of Lemma 5. If $x$ satisfies the condition
$\left({\ref{hypo_r2p}}\right)$ then

$1/\epsilon\leq q:=\frac{x}{2^{\ell}\nu\epsilon}\leq 2/\epsilon.$

Therefore,

$1/\epsilon\leq\left\lceil{{q}}\right\rceil\leq 2/\epsilon.$

The number

$\overline{x}:=2^{\ell}\nu\epsilon\left\lceil{{q}}\right\rceil$

belongs to ${\cal{F}}$ because if $\left\lceil{{q}}\right\rceil<2/\epsilon$
then it fits on definition $\left({\ref{normal}}\right)$ with $e=\ell$ and
$m=\left\lceil{{q}}\right\rceil$, and if
$\left\lceil{{q}}\right\rceil=2/\epsilon$ then

$\overline{x}=2^{\ell+1}\nu\epsilon/\epsilon$

and definition $\left({\ref{normal}}\right)$ holds with $e=\ell+1$ and
$m=1/\epsilon$. We have that

$2^{\ell}\nu\epsilon q=x\leq\overline{x},$

and the definition of ${\mathrm{up}}\\!\left(x\right)$ in
$\left({\ref{up}}\right)$ implies that
${\mathrm{up}}\\!\left(x\right)\leq\overline{x}$. To prove that
${\mathrm{up}}\\!\left(x\right)=\overline{x}$, we now show that if
$y\in{\cal{F}}$ is such that $y\geq x$ then $y\geq\overline{x}$. In fact, if
$y\geq x$ then $y\geq 2^{\ell}\nu$ and this implies that $y$ is normal, that
is,

$y=2^{e_{y}}\nu\epsilon m_{y}\hskip 14.22636pt\mathrm{with}\hskip
14.22636pt1/\epsilon\leq m_{y}<2/\epsilon\hskip 14.22636pt\mathrm{and}\hskip
14.22636pt2^{e_{y}}\nu\epsilon m_{y}\geq x=2^{\ell}\nu\epsilon q.$

This leads to

$2^{e_{y}}m_{y}\geq 2^{\ell}q\ \ \Rightarrow\ \ 2^{e_{y}-\ell}\geq
q/m_{y}>\left({1/\epsilon}\right)/(2/\epsilon)=1/2\Rightarrow e_{y}\geq\ell,$

and $2^{e_{y}-\ell}m_{y}$ is integer. As a result

$2^{e_{y}}m_{y}\geq 2^{\ell}q\ \Rightarrow\ 2^{e_{y}-\ell}m_{y}\geq q\
\Rightarrow\ 2^{e_{y}-\ell}m_{y}\geq\left\lceil{{q}}\right\rceil\ \Rightarrow\
y=2^{e_{y}}\nu\epsilon m_{y}\geq
2^{\ell}\nu\epsilon\left\lceil{{q}}\right\rceil=\overline{x},$

and we are done. $\blacksquare$

## References

  * [1] De Berg. M., van Kreveld, M., Overmans, M., Schwarzkopf,O., Computational Geometry, algorithms and applications, Springer (2008).
  * [2] Dekker, T.J., A floating-point technique for extending the available precision. Numer. Math. 18(3), 224 (1971).
  * [3] Demmel, J., Hida, Y., Fast and accurate floating point summation with application to computational geometry. Numer. Algorithms 37(1–4), 101–112 (2004).
  * [4] Espelid, T. O., On floating point summation. SIAM Review 37, 603–607 (1995).
  * [5] Graillat, S., Louvet, N., Applications of fast and accurate summation in computational geometry. Technical report, Laboratoire LP2A, University of Perpignan, Perpignan, France (2006).
  * [6] Higham, N. J., The accuracy of floating point summation. SIAM J. Sci. Computation 14, 783–799 (1993).
  * [7] Kahan, W. Further remarks on reducing truncation errors. Commun. ACM 8(1), 40 (1965).
  * [8] Lange, M., Oishi, S., A note on Dekker’s FastTwoSum algorithm. Numerische Mathematik, 145(2), 383–403 (2020).
  * [9] Mascarenhas, W.F., Moore: Interval Arithmetic in C++20, In: Barreto G., Coelho R. (eds) Fuzzy Information Processing. NAFIPS 2018. Communications in Computer and Information Science, vol 831, pp 519–529 (2018).
  * [10] Mascarenhas, W.F., Floating point numbers are real numbers, arXiv:1605.09202v2 [math.NA] (2017).
  * [11] Ogita, T., Rump, S.M., Oishi, S, Accurate sum and dot product. SIAM J. Sci. Comput. 26(6), 1955–1988 (2005).
  * [12] Ozaki K., Ogita, T., Oishi S., A robust algorithm for geometric predicate by error-free determinant transformation, Information and Computation 216 3–13 (2012).
  * [13] Priest, D. M., On properties of floating point arithmetics: numerical stability and the cost of accurate computations. Ph.D. Thesis, University of California, Berkeley (1992).
  * [14] Ratschek,H. and Rokne, J., Geometric Computations with Interval and New Robust Methods, Applications in Computer Graphics, GIS and Computational Geometry. Horwood Publishing Chichester (2003).
  * [15] Rump, S., Error estimation of floating-point summation and dot product. BIT Numerical Mathematics 52(1) 201–220 (2012).
  * [16] Shewchuk, J., Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete & Computational Geometry 18:305–363 (1997).
  * [17] Zhu, Y.K., Hayes, W.B., Algorithm 908. ACM Trans. Math. Softw. 37(3), 1–13 (2010).