# Root finding with interval arithmetic

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

###### Abstract

We consider the solution of nonlinear equations in one real variable, the
problem usually called by root finding. Although this is an old problem, we
believe that some aspects of its solution using interval arithmetic are not
well understood, and we present our views on this subject. We argue that
problems with just one variable are much simpler than problems with more
variables, and we should use specific methods for them. We provide an
implementation of our ideas in C++, and make this code available under the
Mozilla Public License 2.0.

##  1 Introduction

Many books and articles explain how to solve the nonlinear equation
${f}\\!\left(x\right)=0$ for $x\in\mathds{R}{}$, and some of them consider the
verified solution of such equation, that is, finding solutions with rigorous
bounds on them. Here we discuss the computation of verified solutions using
interval arithmetic. This is an old subject, but the interval arithmetic
literature is mostly concerned with the solution of the general equation
${f}\\!\left(x\right)=0$ for $x\in{\mathds{R}}^{n}$ using variations of
Newton’s method, since its very beginning [2, 10, 11]. For instance, the third
chapter of [13] gives a nice description of interval Newton’s method and [2,
3] present interesting ways to improve it. Here we focus on the simplest case
$n=1$, because it is very important in practice and there are techniques which
are applicable only for $n=1$. As far as we know, the most detailed discussion
of this particular case is presented in chapter 9 of [4], but we believe that
there is more to be said about this subject than what is presented there.

In this article, $f$ is a function from $\mathds{R}{}$ to $\mathds{R}{}$. We
denote $f$’s $\ell$th derivative by $f^{\left({\ell}\right)}$ and when we
mention such derivatives we assume that they exist. We consider the class

${\mathds{I}}:=\left\\{{\,\mathbb{x}=[\underline{x},\overline{x}],\
\mathrm{with}\
-\infty\leq\underline{x}\leq\overline{x}\leq+\infty}\right\\}\cup\left\\{{\emptyset}\right\\}$

of all closed intervals, including the empty and unbounded ones. For
${\cal{S}}\subset\mathds{R}{}$, we write

${f}\\!\left({\cal{S}}\right):=\left\\{{{f}\\!\left(x\right)\ \mathrm{for}\
x\in{\cal{S}}}\right\\}.$

We assume that we have extensions of $f$ and its derivatives, in the sense
that:

###### Definition 1

We say that a function $\mathbb{f}:{\mathds{I}}\rightarrow{\mathds{I}}$ is an
extension of a function $f:\mathds{R}{}\rightarrow\mathds{R}{}$ in an interval
$\mathbb{x}$ if for every $\mathbb{x}^{\prime}\subset\mathbb{x}$ we have that
${f}\\!\left(\mathbb{x}^{\prime}\right)\subset{\mathbb{f}}\\!\left(\mathbb{x}^{\prime}\right)$.
$\blacktriangle$

We usually identify the point $t\in\mathds{R}{}$ with the interval $[t,t]$,
and for functions $\mathbb{f}:{\mathds{I}}\rightarrow{\mathds{I}}$ we write
${\mathbb{f}}\\!\left(t\right):={\mathbb{f}}\\!\left([t,t]\right)$. The set of
roots is denoted by ${\cal{R}}$, and $r$ is a generic root.

When $f$ is differentiable, the main tool for root finding is Newton’s method

$x_{k+1}=x_{k}-{f}\\!\left(x_{k}\right)/{f}^{\prime}\\!\left(x_{k}\right).$
(1)

This method has a natural extension to interval arithmetic:

$\mathbb{x}_{k+1}=\mathbb{x}_{k}\cap\left({\
t_{k}-{\mathbb{f}}\\!\left(t_{k}\right)/\,{\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)\
}\right),$ (2)

where $t_{k}$ is a point in $\mathbb{x}_{k}$ and $f$ and $f^{\prime}$ are
extensions of $\mathbb{f}$ and $\mathbb{f}^{\prime}$. Traditionally, we
compute $t_{k}$ by rounding the midpoint of $\mathbb{x}_{k}$. The first
question an alert numerical analyst would ask about Equation
$\left({\ref{newton_i}}\right)$ is:

$\mathrm{What\ should\ we\ do\ when\
}{\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)\ \mathrm{contains}\
0\mathrm{?}$ (3)

However, we see few people asking the following question:

$\mathrm{What\ should\ we\ do\ when\
}{\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)\ \mathrm{{\bf does\ not}\
contain}\ 0\mathrm{?}$ (4)

Both questions are quite relevant, but before we answer them we must say that
the interval version of Newton’s method may be implemented without Equation
$\left({\ref{newton_i}}\right)$. Instead, we could write $f$ in the centered
form

${f}\\!\left(t\right)\in{\mathbb{f}}\\!\left(t_{k}\right)+{\mathbb{s}_{t_{k}}}\\!\left(\mathbb{x}_{k}\right)\left({\mathbb{x}_{k}-t_{k}}\right)\hskip
14.22636pt\mathrm{for}\hskip 14.22636ptt\in\mathbb{x}_{k},$ (5)

and replace the extended derivative
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$ by an extension
${\mathbb{s}_{t_{k}}}\\!\left(\mathbb{x}_{k}\right)$ of the slope $s_{t_{k}}$
of $f$ at $t_{k}$. This leads to an improved version of Newton’s method given
by

$\mathbb{x}_{k+1}=\mathbb{x}_{k}\cap\left({\
t_{k}-{\mathbb{f}}\\!\left(t_{k}\right)/\,{\mathbb{s}_{t_{k}}}\\!\left(\mathbb{x_{k}}\right)\
}\right).$ (6)

The slope is defined as

${s_{c}}\\!\left(t\right):={s_{f,c}}\\!\left(t\right):=\left\\{\begin{array}[]{ccc}{f}^{\prime}\\!\left(c\right)&\mathrm{if}&t=c,\\\\[2.84544pt]
\frac{{f}\\!\left(t\right)-{f}\\!\left(c\right)}{t-c}&\mathrm{if}&t\neq
c,\end{array}\right.$ (7)

and centered forms are already mentioned in Moore’s book [10]. They have been
discussed in detail in the interval arithmetic literature [7, 12] and have
generalizations called Taylor forms or Taylor models [9]. Our algorithm uses
only the plain centered form $\left({\ref{center}}\right)$, and in situations
in which practical implementations of such generalizations are available we
would use them only as a tool to obtain more accurate centered forms.

The Mean Value Theorem shows that

$t\leq
c\Rightarrow{s_{c}}\\!\left(t\right)\in{f}^{\prime}\\!\left([t,c]\right)\hskip
28.45274pt\mathrm{and}\hskip 28.45274ptt\geq
c\Rightarrow{s_{c}}\\!\left(t\right)\in{f}^{\prime}\\!\left([c,t]\right)$

and this implies that any extension of $\mathbb{f}^{\prime}$ is an extension
of $s_{t_{k}}$, but there may be better ones, specially when the interval
$\mathbb{x}_{k}$ is not small. For instance, if ${f}\\!\left(t\right):=t^{2}$
then

${s_{0}}\\!\left(t\right)=x\hskip 28.45274pt\mathrm{and}\hskip
28.45274pt{f}^{\prime}\\!\left(t\right)=2t$

and any extension of $f^{\prime}$ yields intervals twice as large than
necessary for an slope.

In practice the gain by replacing derivatives by slopes may not be great,
because usually ${f}^{\prime}\\!\left(r\right)\neq 0$ at the roots $r$ and in
this case, for $t_{k}$ close to $r$, the difference between $s_{t_{k}}$ and
${f}^{\prime}\\!\left(t\right)$ for $t$ in a short interval $\mathbb{x}_{k}$
is not large. Moreover, extensions of derivatives have an important feature
that extensions of slopes do not have: derivative extensions can detect
monotonicity. In other words, if
${\mathbb{f}^{\prime}}\\!\left(\mathbb{x}\right)\subset(0,+\infty)$ then $f$
is increasing in $\mathbb{x}$, but we cannot reach the same conclusion from
the fact that ${s_{t_{k}}}\\!\left(\mathbb{x}\right)\subset(0,+\infty)$. The
ease with which we can use information about monotonicity when we have only
one variable is what makes the case $n=1$ so special, and this is the reason
why this article was written to start with.

We can ask questions similar to $\left({\ref{qd0}}\right)$ and
$\left({\ref{qd1}}\right)$ about the modified Newton’s step
$\left({\ref{newton_t1}}\right)$:

$\mathrm{What\ should\ we\ do\ when\
}{\mathbb{s}_{m_{k}}}\\!\left(\mathbb{x}_{k}\right)\ \mathrm{contains}\
0\mathrm{?}$ (8) $\mathrm{What\ should\ we\ do\ when\
}{\mathbb{s}_{m_{k}}}\\!\left(\mathbb{x}_{k}\right)\ \mathrm{does\ not\
contain}\ 0\mathrm{?}$ (9)

and there are at least two more questions we must ask:

$\mathrm{What\ information\ about\ the\ roots\ our\ algorithm\ should\
provide?}$ (10) $\mathrm{When\ should\ we\ stop\ the\ iterations?}$ (11)

After much thought about the questions above, we devised a root finding
algorithm which combines the interval versions of Newton’s method in Equations
$\left({\ref{newton_i}}\right)$ and $\left({\ref{newton_t1}}\right)$, and
modifies them in order to exploit montonicity. Our algorithm tries to strike a
balance between theory and practice, and takes into account the following
practical issues:

  1. 1.

Usually, the evaluation of interval extensions over intervals is much less
accurate than the evaluation of interval extensions at points, that is, the
width of the interval ${\mathbb{f}}\\!\left(t_{k}\right)$ is much smaller than
the width of the interval ${\mathbb{f}}\\!\left(\mathbb{x}_{k}\right)$. The
same applies to the derivative $\mathbb{f}^{\prime}$.

  2. 2.

Usually, the floating point evaluation of
$d_{k}={f}^{\prime}\\!\left(t_{k}\right)$ yields a reasonable estimate this
derivative, at a much lower cost than the evaluation of the extensions
${\mathbb{f}}^{\prime}\\!\left(t_{k}\right)$ or
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$. The only defect of
$d_{k}$ is that it does not come with a guarantee of its accuracy. As a result
we can, and should, use floating point computations in order to obtain
reasonable estimates (which may turn out to be inaccurate sometimes) and
resort to interval computation only when we absolutely need guarantees about
our results. In particular, the computed
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$ may be very wide even
when the floating point $d_{k}={f}^{\prime}\\!\left(t_{k}\right)$ would lead
to a good Newton step $t_{k+1}=t_{k}-\hat{w}_{k}/d_{k}$, where $\hat{w}_{k}$
is the mid point of $\mathbb{w}_{k}={\mathbb{f}}\\!\left(t_{k}\right)$.

  3. 3.

In interval arithmetic, our goal is to find short intervals $\mathbb{r}$ which
may contain roots and to discard intervals guaranteed not to contain roots.

  4. 4.

The simplest way to ensure that an interval
$\mathbb{r}=[\underline{r},\overline{r}]$ contains a root is to prove that
${f}\\!\left(\underline{r}\right)$ and ${f}\\!\left(\overline{r}\right)$ have
opposite signs. We can do that by evaluating $\mathbb{f}$ at the points
$\underline{r}$ and $\overline{r}$, and as we noted above this evaluation
tends do yield sharp results.

  5. 5.

The simplest way to ensure that an interval $\mathbb{r}$ does not contain a
root is to prove that ${f}\\!\left(\underline{r}\right)$ and
${f}\\!\left(\overline{r}\right)$ are different from zero and have the same
sign and ${f}^{\prime}\\!\left(\mathbb{r}\right)$ does not contain $0$. In
practice, this is not as easy as in the previous item because the computed
${\mathbb{f}}^{\prime}\\!\left(\mathbb{r}\right)$ tends to be inflated by the
usual weaknesses of interval arithmetic. However, when we know that $f$ is
monotonic we can dispense with
${\mathbb{f}}^{\prime}\\!\left(\mathbb{r}\right)$ and check only the values of
$f$ at the end points. This simplifies things immensely for monotonic
functions. Actually, it makes little practical sense to compute
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$ once we already know
that it does not contain 0, and we have this information for all nodes below a
first node $\mathbb{x}_{k}$ in the branch and bound tree which is such that
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$ does not contain 0.

  6. 6.

Multiple roots, i.e., roots $r$ such that ${f}^{\prime}\left(r\right)=0$, are
a major problem and should be handled differently from single roots. Usually,
it is hopeless to try to get them with the same accuracy as single roots, and
the cases in which we have some hope are too specific to deserve attention in
a generic software.

Combining the items above, we propose an algorithm which can be outlined as
follows. For each candidate interval we keep seven additional bits of
information:

  * •

The sign $\sigma_{i}\in\left\\{{-1,0,1}\right\\}$ of $f$ at its infimum,

  * •

The sign $\sigma_{s}\in\left\\{{-1,0,1}\right\\}$ of $f$ at its supremum

  * •

The sign $\sigma_{d}\in\left\\{{-1,0,1}\right\\}$ of $f^{\prime}$ in
$\mathbb{x}$

  * •

A point $t$ in its interior, indicating where it should be split in the Newton
step

  * •

A point $\tilde{t}$, which lies near $\mathbb{x}$.

  * •

The value $\tilde{d}={f}^{\prime}\\!\left(\tilde{t}\right)$.

  * •

The expected sign $\sigma_{t}$ of ${f}\\!\left(t\right)$.

Regarding the signs above, $-1$ means definitely negative, $+1$ means
definitely positive and $0$ means that we do not know the corresponding sign.

We then procedure in the usual branch and bound way:

  1. 1.

If the stack is empty then we are done. Otherwise we pop a tuple

$\left({\mathbb{x}_{k},\sigma_{i},\sigma_{s},\sigma_{d},t_{k},\tilde{t},\tilde{d},\sigma_{t}}\right)$

from the stack.

  2. 2.

If the width of $\mathbb{x}$ is below a given tolerance $\tau_{x}$ then we
insert $\left({\mathbb{x},\sigma_{i},\sigma_{s},\sigma_{d}}\right)$ in a list
of possible solutions and go to item 1.

  3. 3.

We compute ${\mathbb{f}}\\!\left(\mathbb{x}\right)$. If it does not contain
$0$, then drop $\mathbb{x}$ and go to item 1. (Tolerances are discussed in
Section 3.)

  4. 4.

We compute the slope
$\mathbb{s}_{k}:={\mathbb{s}_{t_{k}}}\\!\left(\mathbb{x}_{k}\right)$. If it
does not contain $0$ then we compute the derivative
$\mathbb{d}_{k}:={\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$. If
$\mathbb{d}_{k}$ does not contain $0$ then we change to the specific algorithm
to handle strictly monotonic functions described in Section 5, and once we are
done we go to item 1. Otherwise, we replace $\mathbb{s}_{k}$ by
$\mathbb{s}_{k}\cap\mathbb{d}_{k}$ and continue.

  5. 5.

We compute ${\mathbb{f}}\\!\left(t_{k}\right)$ and check whether it contains
zero or is too close to zero. If it does then use the algorithm described in
Section 3 to handle $\mathbb{x}_{k}$ and go to item 1.

  6. 6.

If ${\mathbb{f}}\\!\left(t_{k}\right)$ is not too close to zero then and apply
the version of the interval Newton step described in Section 2, obtaining at
most two intervals $\mathbb{x}_{k+1}$ and $\mathbb{x}_{k+1}^{\prime}$. If
there are no $\mathbb{x}_{k+1}$ then we drop $\mathbb{x}_{k}$ and go to item
1.

  7. 7.

If there is just one $\mathbb{x}_{k+1}$ then we check whether the sign
$\sigma_{t}$ matches the sign of ${\mathbb{f}}\\!\left(\mathbb{x}_{k}\right)$.
If it does not then we set $\sigma_{t}$ to zero, take $t_{k+1}$ as the mid
point of $\mathbb{x}_{k+1}$, obtain the extra information for
$\mathbb{x}_{x+1}$ and push it on the stack and go to item 1. Otherwise, we
compute $d={f}^{\prime}\\!\left(t_{k}\right)$ using floating point arithmetic
and use $\tilde{d}$ and $\tilde{t}$ to check whether the corrected Newton step
$t_{k+1}$ described in Section 4 lies in $\mathbb{x}_{k+1}$. If it does then
we use this $t_{k+1}$, otherwise we take $t_{k+1}$ as the midpoint of
$\mathbb{x}_{k+1}$. We then obtain the additional information for
$\mathbb{x}_{k+1}$, push it on the stack and go to item 1.

  8. 8.

If there are two $\mathbb{x}_{k+1}$s, then we check whether the width of
$\mathbb{x}_{k}$ lies below the cluster threshold $\tau_{c}$. If it does then
we deem $\mathbb{x}_{k}$ to contain a cluster of zeros, and proceed as in item
2. Otherwise we apply the same procedure as in item 7 to $\mathbb{x}_{k+1}$
and $\mathbb{x}_{k+1}^{\prime}$ and go to item 1.

We implemented the algorithm outlined above in C++, using our Moore interval
arithmetic library [8]. This code is available with the arxiv version of this
article. It is distributed under the Mozilla Public License 2.0.
Unfortunately, the code is much more involved than the summary above, because
it must deal with many technical details which were not mentioned in this
article in order not to make it longer and more complex than it already is.

In the rest of the article we discuss in more detail several aspects of the
algorithm outline above. We start with Section 2, in which we present a
version of Newton’s method for interval arithmetic. This version is similar to
the ones found in the literature, but it is slightly different because it
ensures that $f$ is different from zero at some of the extreme points of the
new intervals, and computes the signs of $f$ at these extremes at a low extra
cost. As a result, our algorithm yields not only intervals containing the
roots but also the sign of $f$ at the extreme of such intervals, and these
signs certify the existence of a root in the interval when they are different
from each other. In Section 3 we make some comments about interval arithmetic
in general and discuss the thorny subject of tolerances, which are unavoidable
for defining stopping criterion for interval root solvers. Section 4 presents
yet another version of the classic Newton’s method for exact, “point
arithmetic.” This “point version” is the motivation for the interval version
of Newton’s method for monotone functions presented in Section 5. Finally, in
Section 6 we discuss the important subject of testing.

##  2 The Interval version of Newton’s step

This section is about the interval version of Newton’s step

$\mathbb{x}_{k+1}=\mathbb{x}_{k}\cap\left({\
t_{k}-{\mathbb{f}}\\!\left(t_{k}\right)/\,\mathbb{d}_{k}\ }\right),$ (12)

where $\mathbb{d}_{k}=[\underline{d}_{k},\overline{d}_{k}]$ can be either the
derivative ${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$ or the slope
${\mathbb{s}_{t_{k}}}\\!\left(\mathbb{x_{k}}\right)$. Here we make the
simplifying assumption that the interval

$\mathbb{w}_{k}:={\mathbb{f}}\\!\left(t_{k}\right):=[\underline{w}_{k},\
\overline{w}_{k}]$

does not contain $0$, answering to the questions $\left({\ref{qd0}}\right)$,
$\left({\ref{qd1}}\right)$, $\left({\ref{qs0}}\right)$ and
$\left({\ref{qs1}}\right)$ in the introduction in this case. By replacing $f$
by $-f$ if necessary, we can assume that $\overline{w}_{k}<0$, and we make
this simplifying assumption from now on.

The answer to questions $\left({\ref{qd0}}\right)$ and
$\left({\ref{qs0}}\right)$, in which case $\mathbb{d}_{k}$ contains $0$, is
illustrated in Figure 2. In this Figure, the inclined lines have equations

$w=\overline{w}_{k}+\underline{d}_{k}\left({t-t_{k}}\right)\hskip
28.45274pt\mathrm{and}\hskip
28.45274ptw=\overline{w}_{k}+\overline{d}_{k}\left({t-t_{k}}\right),$ (13)

and by intersecting these lines with the axis $w=0$ we obtain a better bound
on ${\cal{R}}\cap\mathbb{x}_{k}$ (recall that ${\cal{R}}$ is the set of
roots.) There are three possibilities for our bounds on
${\cal{R}}\cap\mathbb{x}_{k}$ after we take in to account the intersections
above. We may find that:

  * •

${\cal{R}}\cap\mathbb{x}_{k}=\emptyset$. In this case we drop
$\mathbb{x}_{k}$.

  * •

${\cal{R}}\cap\mathbb{x}_{k}$ is contained in a single interval
$\mathbb{r}=[\underline{r},\overline{r}]$, as in the second and third cases in
Figure 1. In this case we take $\mathbb{x}_{k+1}=\mathbb{r}$.

  * •

${\cal{R}}\cap\mathbb{x}_{k}$ is contained in the union of two intervals
$\mathbb{r}_{1}$ and $\mathbb{r}_{2}$ such that
$\overline{r}_{1}\leq\underline{r}_{2}$, as in the last case of Figure 1.

${\cal{R}}\cap\mathbb{x}_{k}$ is
empty$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$\underline{r}$$\overline{r}$$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$\underline{r}$$\overline{r}$$\overline{w}_{k}$$\underline{w}_{k}$$t_{k}$$\underline{r}_{1}$$\overline{r}_{1}$$\underline{r}_{2}$$\overline{r}_{2}$
Figure 1: The case $0\in\mathbb{d}$, $\overline{w}<0$, with
$-\underline{d}_{k}$ and $\overline{d}_{k}$ large and small.

There is nothing new in our approach up to this point. We just explained with
a picture what most people do algebraically. Our slight improvement comes in
the way we compute the intersections in Equations
$\left({\ref{inter1}}\right)$. These intersections are

$\overline{r}_{1}:=t_{k}-\overline{w}_{k}/\underline{d}_{k}\hskip
28.45274pt\mathrm{and}\hskip
28.45274pt\underline{r}_{2}:=t_{k}-\overline{w}_{k}/\overline{d}_{k},$

and, as usual in the Moore Library [8], we propose that $\overline{r}_{1}$ and
$\underline{r}_{2}$ be computed with the rounding mode upwards, using the
following expressions:

$\displaystyle\overline{r}_{1}$ $\displaystyle:=$ $\displaystyle
t_{k}+\overline{q}_{k}\hskip 14.22636pt\mathrm{for}\hskip
8.5359pt\overline{q}_{k}:=\overline{w}_{k}/\left({-\underline{d}_{k}}\right),$
(14) $\displaystyle\underline{r}_{2}$ $\displaystyle:=$ $\displaystyle-
u_{k}\hskip 28.45274pt\mathrm{for}\hskip
8.5359ptu_{k}:=\underline{q}_{k}-t_{k}\hskip 8.5359pt\mathrm{and}\hskip
8.5359pt\underline{q}_{k}:=\overline{w}_{k}/\overline{d}_{k}.$ (15)

This order of evaluation ensure that $\overline{r}_{1}$ and
$\underline{r}_{2}$ are rounded in the correct direction, so that we do not
risk loosing roots. The expressions in Equation
$\left({\ref{computedr1}}\right)$ and $\left({\ref{computedr2}}\right)$ are so
simple that we find whether $\overline{r}_{1}$ and $\underline{r}_{2}$ where
computed exactly by checking whether

$\overline{r}_{1}-t_{k}=\overline{q}_{k}\hskip 14.22636pt\mathrm{and}\hskip
14.22636pt\overline{q}_{k}\underline{d}_{k}=\overline{w}_{k}$ (16)

and

$u_{k}+t_{k}=\underline{q}_{k}\hskip 14.22636pt\mathrm{and}\hskip
14.22636pt\underline{q}_{k}\overline{d}_{k}=\overline{w}_{k}.$ (17)

If we find that $\overline{r}_{1}$ was computed exactly then we increase it to
the next floating point number. If $\underline{r}_{2}$ was computed exactly
then we decrease it to the previous floating point number. By doing so, we
ensure that ${f}\\!\left(\overline{r}_{1}\right)<0$ and
${f}\\!\left(\underline{r}_{2}\right)<0$, without computing
${\mathbb{f}}\\!\left(\overline{r}_{1}\right)$ or
${\mathbb{f}}\\!\left(\underline{r}_{2}\right)$. We should also mention that
it is possible to prove that, even after the rounding, incrementing and
decrementing steps above, $\overline{r}_{1}\leq t_{k}\leq\underline{r}_{2}$
and there is no risk of $\overline{r}_{1}$ crossing $\underline{r}_{2}$.

Regarding the cost of all this, note that in the Moore library the rounding
mode is always upwards. Therefore, there is no cost associated with changing
rounding modes. Moreover, the expensive operations in Equations
$\left({\ref{computedr1}}\right)$ and $\left({\ref{computedr2}}\right)$ are
the divisions, and the extra couple of sums and multiplications do not
increase the overall cost of computing the intersections by much. In fact, if
we take into account branch prediction and speculative execution and the fact
that computations are usually inexact, the extra cost due to the verification
in Equations $\left({\ref{veri1}}\right)$ and $\left({\ref{veri2}}\right)$ and
the decrement of $\overline{r}_{1}$ and the increment of $\underline{r}_{2}$
is likely to be minimal. We would not be surprised if by performing the
verification in Equations $\left({\ref{veri1}}\right)$ and
$\left({\ref{veri2}}\right)$ above the code would be faster than blindly
incrementing $\overline{r}_{1}$ and decrementing $\underline{r}_{2}$ (we did
not have the time to check this.)

Finally, by using symmetry, we can reduce the analysis of questions
$\left({\ref{qd1}}\right)$ and $\left({\ref{qs1}}\right)$ in the case in which
${\mathbb{f}}\\!\left(t_{k}\right)$ does not contain $0$ to the cases
described in Figure 2. In this Figure, the inclined lines have equations

$w=\overline{w}_{k}+\overline{d}_{k}\left({t-t_{k}}\right)\hskip
28.45274pt\mathrm{and}\hskip
28.45274ptw=\underline{w}_{k}+\underline{d}_{k}\left({t-t_{k}}\right).$

and we may either have no intersection or one intersection. We can use the
same technique described above for find intersections which are correctly
rounded and such that ${f}\\!\left(u\right)$ is different from zero in the new
extreme points $u$, without evaluating ${\mathbb{f}}\\!\left(u\right)$ (the
old extreme points stay as they were.)

${\cal{R}}\cap\mathbb{x}_{k}$ is
empty$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$\underline{r}$$\overline{r}$$t_{k}$$\overline{w}_{k}$$\underline{w}_{k}$$\underline{r}$$\overline{r}$
Figure 2: The case $\underline{d}_{k}>0$ and $\overline{w}_{k}<0$, with
$\underline{d}_{k}$ and $\overline{d}_{k}$ large and small.

##  3 What we can expect from Interval Arithmetic

In order to appreciate the content of this article, one must realize that
interval arithmetic is quite different from floating point arithmetic, and it
is used for different purposes. Floating point methods are usually faster than
interval methods, but do not provide rigorous bounds on their results. We use
interval arithmetic when we want global solutions to our problems, with
rigorous bounds on them, and we are willing to pay more for that.

$\underline{r}$$\overline{r}$$\underline{r}^{\prime}$$\overline{r}^{\prime}$${f}\\!\left(t\right)$${\underline{\mathbb{f}}}\\!\left(t\right)$${\overline{\mathbb{f}}}\\!\left(t\right)$$r$
Figure 3: If $\mathbb{f}$ is all we know, then the best estimate we can hope
for $r$ is the interval $\mathbb{r}=[\underline{r},\overline{r}]$. Knowing
that ${\underline{\mathbb{f}}}^{\prime}\\!\left(t\right)>0$ for all
$t\in\mathbb{r}$ we can improve this to
$\mathbb{r}^{\prime}=[\underline{r}^{\prime},\overline{r}^{\prime}]$. This is
the best we can do if $\mathbb{f}$ an $\mathbb{f}^{\prime}$ are all the
information we have about $f$.

Interval arithmetic is fundamentally limited. When we evaluate the extension
$\mathbb{f}$ at a point $t$ obtain an interval
$\mathbb{w}={\mathbb{f}}\\!\left(t\right)$ containing ${f}\\!\left(t\right)$.
The width of $\mathbb{w}$ depends on the quality of the implementation of
$\mathbb{f}$ and the underlying arithmetic. When $f$ is not too complex and
${f}\\!\left(t\right)$ is ${O}\\!\left(1\right)$, we expect that the width of
$\mathbb{w}$ to be ${O}\\!\left(\epsilon\right)$, where $\epsilon$ is the
machine precision. It is also reasonable to expect that the functions
$\underline{\mathbb{f}},\overline{\mathbb{f}}:\mathbb{x}\rightarrow\mathds{R}{}$
will be as in Figure 3 when $r\in\mathbb{x}$ is a root of $f$.

When the evaluation of $f$ using interval arithmetic leads to wide intervals,
we can (and do) use centered forms and Taylor models to reduce the width of
${\mathbb{f}}\\!\left(\mathbb{x}\right)$, but this will not free us from
facing the reality that the evaluation of interval extensions is usually
imperfect. This is a fact, and we need to be prepared to handle not only the
monotone situation described in Figure 3, but also the noisy situation
described in Figure 4, which often occurs near a double root, and gets worse
near roots of higher multiplicities.

${f_{1}}\\!\left(t\right)$${f_{2}}\\!\left(t\right)$${\underline{\mathbb{f}}}\\!\left(t\right)$${\overline{\mathbb{f}}}\\!\left(t\right)$$\tau_{w}$$\underline{x}$$\underline{c}$$\tau_{c}$$z$$\tau_{c}$$\tau_{c}$$\overline{c}$$\overline{x}$
Figure 4: The hardest case our algorithm faces: a cluster of zeros below the
accuracy of $\mathbb{f}$ and $\mathbb{f}^{\prime}$. There are no silver
bullets in such cases, and we must find reasonable ways to cope with the
$\mathbb{f}$ and $\mathbb{f}^{\prime}$ which are provided to us, not the ones
that we dream about. Using the tolerances $\tau_{c}$ and $\tau_{w}$ described
in the text is a reasonable compromise to handle such situations (which is far
from perfect.)

Another fundamental fact about interval arithmetic is that it is a pessimistic
theory: it always considers the worst scenario, and it would be inconsistent
if it were otherwise. As consequence, of this fundamental fact, the functions
$f_{1}$ and $f_{2}$ are indistinguishable by their common extension
$\mathbb{f}$, that is, if $\mathbb{f}$ is all the information that we have,
then all our conclusions about $f_{1}$ must also apply to $f_{2}$. This forces
us to be very conservative, and include all possible roots for $f_{1}$ as
candidates for roots for $f_{2}$ too, and vice versa. We emphasize this point
because we found that it is not always clear in the minds of people trying to
implement roots finders like the one we propose here, and this makes them
under estimate the effort required for this task in real life.

Taking into account all that was said since the beginning of this section, we
reached the conclusion that we should use three tolerances to try to handle
the problems caused by the very nature of interval arithmetic. First, we
should have a tolerance $\tau_{x}$ to handle “normal” roots $r$ like the one
in Figure 3, at which ${f}^{\prime}\\!\left(r\right)\neq 0$, so that we
consider intervals around such roots with width below $\tau_{x}$ to be good
enough. Clusters of zeros of $\mathbb{f}$, as in Figure 4, are a different
kind of beast, as the experience with other problems in interval arithmetic
has show [5, 6]. By using the same tolerance $\tau_{x}$ for them we can end up
with literally thousands of small intervals as candidates for containing
roots, in the favorable case in which the program actually ends and returns
something. We say this from our experience with our own interval root finders
as well as with interval root finders developed by other people.

Our algorithm asks the user to also provide a “cluster tolerance” $\tau_{c}$,
which could be something of the order of $\sqrt{\tau_{x}}$, and a tolerance
$\tau_{w}$ so that function values $w={f}\\!\left(t\right)$ with magnitude
below $\tau_{w}$ are deemed to be negligible. We take the liberty of
increasing $\tau_{w}$ if we find it to be smaller than $16$ times the maximum
width of the intervals ${\mathbb{f}}\\!\left(t\right)$ which compute along the
execution of algorithm. Contrary to what is proposed in [4], we believe that
all tolerances should be absolute, and not relative (and users are free to
choose which approach suits their needs best.) Using these tolerances, we can
implement Algorithm 1 which “expands” a zero $z\in\mathbb{x}$. By applying
this algorithms to the point $z$ and the interval
$\mathbb{x}=[\underline{x},\overline{x}]$ in Figure 4 we would identify the
cluster $[\underline{c},\overline{c}]$ and submit the intervals
$[\underline{x},\underline{c}]$ and $[\overline{c},\overline{x}]$ to further
examination. Finally, we must say that the actual algorithms for zero
expansion in the C++ code are more involved than Algorithm 1, but the details
are too cumbersome to be presented here.

Algorithm 1 Zero expansion

procedure expand_zero($z$, $\mathbb{x}$)

$\overline{c}\leftarrow z$

while $\overline{c}+\tau_{c}\leq\overline{x}$ do

$\mathbb{w}\leftarrow{\mathbb{f}}\\!\left(\overline{c}+\tau_{c}\right)$

$\tau_{w}\leftarrow\max\left\\{{16\
{\mathrm{wid}}\\!\left(\mathbb{w}\right),\tau_{w}}\right\\}$

if $\underline{w}\geq\tau_{w}$ or $\overline{w}\leq-\tau_{w}$ then

break

$\overline{c}\leftarrow\overline{c}+\tau_{c}$

$\underline{c}\leftarrow z$

while $\underline{c}-\tau_{c}\geq\underline{x}$ do

$\mathbb{w}\leftarrow{\mathbb{f}}\\!\left(\underline{c}-\tau_{c}\right)$

$\tau_{w}\leftarrow\max\left\\{{16\
{\mathrm{wid}}\\!\left(\mathbb{w}\right),\tau_{w}}\right\\}$

if $\underline{w}\geq\tau_{w}$ or $\overline{w}\leq-\tau_{w}$ then

break

$\underline{c}\leftarrow\underline{c}-\tau_{c}$

if $\underline{c}>\underline{x}$ then

if $\overline{c}<\overline{x}$ then

return $[\underline{x},\underline{c}],\ \ [\underline{c},\overline{c}],\ \
[\overline{c},\overline{x}]$

else

return $[\underline{x},\underline{c}],\ \ [\underline{c},\overline{c}],\ \
\emptyset$

else

if $\overline{c}<\overline{x}$ then

return $\emptyset,\ \ [\underline{c},\overline{c}],\ \
[\overline{c},\overline{x}]$

else

return $\emptyset,\ \ [\underline{c},\overline{c}],\ \ \emptyset$

end procedure

##  4 The modified Newton’s method

As mentioned in the introduction, the main tool for root finding is Newton’s
method

$t_{k+1}=t_{k}-{f}\\!\left(t_{k}\right)/{f}^{\prime}\\!\left(t_{k}\right).$
(18)

This method was devised to find a “point approximation” to the root $r$, that
is, we for $\tilde{r}$ such that ${\left|{r-\tilde{r}}\right|}$ is small. It
is our opinion that the goal of interval arithmetic is a bit different: we
look for a short interval $\mathbb{r}$ such that $r\in\mathbb{r}$. Of course,
in the end both goals amount to the same, but they suggest slightly different
perspectives. We can rephrase the interval arithmetic goal as finding a short
interval $\mathbb{r}=[\underline{r},\overline{r}]$ such that
${f}\\!\left(\underline{r}\right)$ and ${f}\\!\left(\overline{r}\right)$ have
opposite signs (modulo degenerate cases in which
${f}\\!\left(\underline{r}\right)=0$ or ${f}\\!\left(\overline{r}\right)=0$.)
From this perspective, we believe that we should modify the classic Newton
step $\left({\ref{newton2}}\right)$ so that it produces iterates $t_{k}$ such
that the signs of $w_{k}:={f}\\!\left(t_{k}\right)$ alternate, in order to
ensure that the interval

$\mathbb{r}_{k}:=[\underline{r}_{k},\overline{r}_{k}]\hskip
14.22636pt\mathrm{for}\hskip
14.22636pt\underline{r}_{k}:=\min\left\\{{t_{k},t_{k+1}}\right\\}\hskip
14.22636pt\mathrm{and}\hskip
14.22636pt\overline{r}_{k}:=\max\left\\{{t_{k},t_{k+1}}\right\\}$

always contain a root.

$w_{k}$$t_{k}$$n_{k+1}$$w={f}\\!\left(t\right)$$t_{k+1}:=n_{k}-h_{k}s_{k}^{2}/d_{k}$$n_{k+1}:=t_{k}+s_{k}$$s_{k}:=w_{k}/d_{k}$$d_{k}:={f}^{\prime}\\!\left(t_{k}\right)$$h_{k}\approx{f}^{\prime\prime}\\!\left(t_{k}\right)<0$$w_{k+1}$$w={f}\\!\left(t\right)$$w_{k}$$t_{k}$$t_{k+1}=n_{k}+1$$w_{k+1}$
Figure 5: The modified Newton step when $w_{k}={f}\\!\left(t_{k}\right)<0$ and
$d_{k}={f}^{\prime}\\!\left(t_{k}\right)>0$. At left, $f$ is concave and we
modify the Newton iterate $n_{k+1}$ in order to ensure that $w_{k+1}$ and
$w_{k}$ have opposite signs. At right, $f$ is convex and no modification is
needed.

A simple modification with this purpose is described in Figure 5. This Figure
illustrates only the case in which $w_{k}<0$ and $d_{k}>0$, but cases in which
$w_{k}>0$ or $d_{k}<0$ are similar and can be reduced to the case $w_{k}<0$
and $d_{k}>0$ by replacing $f$ by $-f$ or $t$ by $-t$. A simple way to obtain
a reasonable approximation $h_{k}$ for the second derivative mentioned in
Figure 5 is to use the quotient

$h_{k}=\frac{d_{k-m}-d_{k}}{t_{k-m}-t_{k}},$

with data $d_{k-m}$ and $t_{k-m}$ from previous iterations.

Motivated by Figure 5, we propose Algorithm 2 below for finding a short
interval $\mathbb{r}$ containing the single root of $f$ in the interval
$\mathbb{x}$ using exact arithmetic, in the case in which
${f}\\!\left(\underline{x}\right)<0<{f}\\!\left(\overline{x}\right)$ and
${f}^{\prime}\\!\left(t\right)>0$ for $t\in\mathbb{x}$. It is not difficult to
prove that Algorithm 2 has the same properties as the many other versions of
Newton’s method that we find in the literature:

  * •

When the model

${f}\\!\left(t\right)\approx{f}\\!\left(t_{k}\right)+{f}^{\prime}\\!\left(t_{k}\right)\left({t-t_{k}}\right)+{f}^{\prime\prime}\\!\left(t\right)\left({t-t_{k}}\right)^{2}/2$
(19)

is inaccurate for $t\in\mathbb{x}$, we fallback into a bisection step. This
guarantees that eventually the interval $\mathbb{x}$ will become very short,
and the analysis in the next item will apply

  * •

When the quadratic model $\left({\ref{quadratic}}\right)$ is accurate in
$\mathbb{x}$, and the interval $\mathbb{x}$ is very short, $h$ will be a good
approximation of ${f}^{\prime\prime}\\!\left(t_{k}\right)$ and our
modification in the Newton step will ensure that bisection will stop being
necessary, the signs of the $w_{k}$ will alternate and the iterates will
converge at the same quadratic rate as the classic method.

Algorithm 2 Find a short interval $\mathbb{r}$ containing the root $r$ of $f$
with ${f}^{\prime}\\!\left(t\right)>0$.

procedure exact_newton($\mathbb{x}$, $f$, $\tau_{x}$, $\tau_{w}$)

$\underline{w}\leftarrow{f}\\!\left(\underline{x}\right),\ \
\overline{w}\leftarrow{f}\\!\left(\overline{x}\right),\ \
\underline{d}\leftarrow\mathrm{nan},\ \ \overline{d}\leftarrow\mathrm{nan},\ \
d\leftarrow\rm{nan},\ \ t_{d}\leftarrow\mathrm{nan},\ \ h\leftarrow 0$

if $\underline{w}>0$ or $\overline{w}<0$ then

return $\emptyset$

regular case:

if $\left({\overline{x}-\underline{x}<\tau_{x}}\right)\ $ or $\
\left({\overline{w}-\underline{w}<\tau_{w}}\right)$ then

return $\mathbb{x}$

if $-\underline{w}\leq\overline{w}$ then

if $!{\mathrm{isnan}}\\!\left(\underline{d}\right)$ then

goto bissection

$t_{k}\leftarrow\underline{x},\ \
d_{k}\leftarrow{f}^{\prime}\\!\left(t_{k}\right)$

if $!{\mathrm{isnan}}\\!\left(d\right)$ and $t_{d}\neq t_{k}$ then

$h\leftarrow\left({d-d_{k}}\right)/\left({t_{d}-t_{k}}\right),\ \ d\leftarrow
d_{k},\ \ t_{d}\leftarrow t_{k}$

$\underline{d}\leftarrow d_{k},\ \ d\leftarrow d_{k},\ \ t_{d}\leftarrow
t_{k},\ \ s_{k}\leftarrow-\underline{w}/d_{k}$

if $h<0$ then

$s_{k}\leftarrow s_{k}-hs_{k}^{2}/d_{k}$

if $2s_{k}>\overline{x}-\underline{x}$ then

goto bissection

$t_{k+1}\leftarrow t_{k}+s_{k},\ \
w_{k+1}\leftarrow{f}\\!\left(t_{k+1}\right)$

if $w_{k+1}\geq 0$ then

$\overline{w}\leftarrow w_{k+1},\ \ \overline{x}\leftarrow t_{k+1},\ \
\overline{d}\leftarrow\mathrm{nan}$

goto regular case

else

$\underline{w}\leftarrow w_{k+1},\ \ \underline{x}\leftarrow t_{k+1},\ \
\underline{d}\leftarrow\mathrm{nan}$

goto bissection

else

the case $-\underline{w}>\overline{w}$ is analogous to
$-\underline{w}\leq\overline{w}$

bissection:

$t_{k+1}\leftarrow\left({\underline{x}+\overline{x}}\right)/2,\ \
w_{k+1}\leftarrow{f}\\!\left(t_{k+1}\right)$

if $w_{k+1}>0$ then

$\overline{w}\leftarrow w_{k+1},\ \ \overline{x}\leftarrow t_{k+1},\ \
\overline{d}\leftarrow\mathrm{nan}$

else

$\underline{w}\leftarrow w_{k+1},\ \ \underline{x}\leftarrow t_{k+1},\ \
\underline{d}\leftarrow\mathrm{nan}$

goto regular case

end procedure

##  5 The Monotonic Method

This Section presents an interval Newton’s method for functions with
${f}^{\prime}\\!\left(t\right)>\kappa>0$ for $t\in\mathbb{x}$ (in order to
used for functions with ${f}^{\prime}\\!\left(t\right)<-\kappa<0$ we can
simply replace $f$ by $-f$. The method is a simple adaptation of the “point
method” presented in Section 4, taking into account the fundamental
differences between exact (point) arithmetic and interval arithmetic. Among
such differences we have the following ones:

  * •

In interval arithmetic, for a strictly increasing function, we may have that
$0\in{\mathbb{f}}\\!\left(t\right)$ for many $t$’s. In exact arithmetic this
can only happen for one value of $t$.

  * •

In interval arithmetic we have a natural measure for the error in the the
values of $f$, given by the maximum width of the intervals
${\mathbb{f}}\\!\left(t_{k}\right)$. This allows us to correct a poor choice
of the tolerance $\tau_{w}$ by the user and dispense with the cluster
tolerance $\tau_{c}$ mentioned in Section 3.

As we mentioned in the introduction, the increasing case is much simpler than
the general one. In this case, we do not need to compute neither
${\mathbb{f}}\\!\left(\mathbb{x}_{k}\right)$ nor
${\mathbb{f}}^{\prime}\\!\left(\mathbb{x}_{k}\right)$. We only need the
interval ${\mathbb{f}}\\!\left(t_{k}\right)$ and the floating point number
${f}^{\prime}\\!\left(t_{k}\right)$, and this entails a considerable reduction
in the cost of the steps. The analysis of degenerate cases is also much
simpler in the increasing case. These ideas are explored in Algorithm 3, which
is presented after the bibliography. For brevity, we omit the expand_zero
label in this code, but it is similar to the zero expansion algorithms
presented in Section 3.

Algorithm 3 Find a short interval $\mathbb{r}$ containing the root $r$ of $f$
with
${\underline{\mathbb{f}}}^{\prime}\\!\left(\mathbb{x}\right)\geq\kappa>0$.

procedure increasing_interval_newton($\mathbb{x}$, $f$, $\tau_{x}$,
$\tau_{w}$)

$\mathbb{w}_{-}\leftarrow{\mathbb{f}}\\!\left(\underline{x}\right),\ \
\mathbb{w}_{+}\leftarrow{\mathbb{f}}\\!\left(\overline{x}\right),\ \
\underline{d}\leftarrow\mathrm{nan},\ \ \overline{d}\leftarrow\mathrm{nan},\ \
d\leftarrow\rm{nan},\ \ t_{d}\leftarrow\mathrm{nan},\ \ h\leftarrow 0$

if $\underline{w}_{-}>0$ or $\overline{w}_{+}<0$ then

return $\emptyset$

$\tau_{w}\leftarrow\max\left\\{{16\
{\mathrm{wid}}\\!\left(\mathbb{w}_{-}\right),16\
{\mathrm{wid}}\\!\left(\mathbb{w}_{+}\right),\tau_{w}}\right\\}$

if $\overline{w}_{-}\geq 0\ \textbf{then}$

$t_{z}\leftarrow\underline{x}$ and goto expand zero

if $\underline{w}_{+}\leq 0\ \textbf{then}$

$t_{z}\leftarrow\overline{x}$ and goto expand zero

regular case:

if $\left({\overline{x}-\underline{x}<\tau_{x}}\right)\ $ or $\
\left({\overline{w}_{+}-\underline{w}_{-}<\tau_{w}}\right)$ then

return $\mathbb{x}$

if $-\underline{w}_{-}\leq\overline{w}_{+}$ then

if $!{\mathrm{isnan}}\\!\left(\underline{d}\right)$ then

goto bissection

$t_{k}\leftarrow\underline{x},\ \
d_{k}\leftarrow\max\left\\{{\kappa,{f}^{\prime}\\!\left(t_{k}\right)}\right\\}$

if $!{\mathrm{isnan}}\\!\left(d\right)$ and $t_{d}\neq t_{k}$ then

$h\leftarrow\left({d-d_{k}}\right)/\left({t_{d}-t_{k}}\right),\ \ d\leftarrow
d_{k},\ \ t_{d}\leftarrow t_{k}$

$\underline{d}\leftarrow d_{k},\ \ d\leftarrow d_{k},\ \ t_{d}\leftarrow
t_{k},\ \ s_{k}\leftarrow-\underline{w}_{-}/d_{k}$

if $h<0$ then

$s_{k}\leftarrow s_{k}-hs_{k}^{2}/d_{k}$

if $2s_{k}>\overline{x}-\underline{x}$ then

goto bissection

$t_{k+1}\leftarrow t_{k}+s_{k},\ \
{\mathbb{w}_{k+1}\leftarrow{\mathbb{f}}\\!\left(t_{k+1}\right),\ \
\tau_{w}\leftarrow\max\left\\{{16\
{\mathrm{wid}}\\!\left(\mathbb{w}_{k+1}\right),\tau_{w}}\right\\}}$

if $\underline{w}_{k+1}>0$ then

$\mathbb{w}_{+}\leftarrow\mathbb{w}_{k+1},\ \ \overline{x}\leftarrow t_{k+1},\
\ \overline{d}\leftarrow\mathrm{nan}$ and goto regular case

else

if $\overline{w}_{k+1}<0$ then

$\underline{w}_{-}\leftarrow\mathbb{w}_{k+1},\ \ \underline{x}\leftarrow
t_{k+1},\ \ \underline{d}\leftarrow\mathrm{nan}$

goto bissection

else

$t_{z}\leftarrow t_{k+1}$ and goto expand zero

else

the case $-\underline{w}_{-}>\overline{w}_{+}$ is analogous to
$-\underline{w}_{-}\leq\overline{w}_{+}$

bissection:

$t_{k+1}\leftarrow\left({\underline{x}+\overline{x}}\right)/2,\ \
\mathbb{w}_{k+1}\leftarrow{\mathbb{f}}\\!\left(t_{k+1}\right)$

if $\underline{w}_{k+1}>0$ then

$\mathbb{w}_{+}\leftarrow\mathbb{w}_{k+1},\ \ \overline{x}\leftarrow t_{k+1},\
\ \overline{d}\leftarrow\mathrm{nan}$

else

if $\overline{w}_{k+1}>0$ then

$\mathbb{w}_{-}\leftarrow\mathbb{w}_{k+1},\ \ \underline{x}\leftarrow
t_{k+1},\ \ \underline{d}\leftarrow\mathrm{nan}$

else

$t_{z}\leftarrow t_{k+1}$ and goto expand zero

goto regular case

end procedure

##  6 Testing

In practice, it is quite difficult to implement the algorithm presented here
due to the need of much attention to detail. It is quite likely that the first
attempts at such an implementation will fail, as ours did. Even worse, it is
also likely that the failure will be undetected, specially if the details
involved in the actual coding of the algorithm are considered to be of little
relevance. As a result, in order to have some assurance of the reliability of
our code in practice it is necessary to have a good suit of tests. Readers
should not underestimate this point (and should not underestimate our remark
that they should not underestimate this point, and, recursively…)

It is also important to realize that a good test suit may reveal problems not
only in the code but also in the theory upon which it is based. A good example
of this point is the choice of the stopping criterion presented in chapter 9
of [4]. When using that stopping criterion for the simple polynomial

${f}\\!\left(x\right)=\left({x-1}\right)\left({x-2}\right)\left({x-3}\right)\left({x-4}\right)\left({x-5}\right)$
(20)

in the interval $\mathbb{x}=[1,5]$, we would obtain that the solution set
${\cal{R}}$ is the whole interval $\mathbb{x}$, regardless of the tolerances
provided by the user, and a good set of tests could detect this problem. The
authors of [4] do mention that their stopping criterion may be problematic in
some rare cases, and that users would detect this kind of problem afterwards.
Readers may think that the polynomial in Equation $\left({\ref{bughw}}\right)$
is one of such cases, but we beg to differ: we believe that a proper algorithm
should get the roots of a function as simple as the one in Equation
$\left({\ref{bughw}}\right)$ with high accuracy without user intervention, and
if it does not then it must be fixed. Asking users to analyze the result may
not be an option for instance when our algorithm is being used thousands of
times as part of a larger routine. In this scenario, users will not have the
chance to look at the result of each call of the algorithm, and they will need
either to rely on the algorithm to provide good answers, or to write their own
code to do what the algorithm should have done for them.

The example in Equation $\left({\ref{bughw}}\right)$ can be though as a mild
version of Wilkinson’s polynomial, and we believe that we can obtain a good
test suit by adding multiple roots to polynomials similar to it. In our
opinion, a robust root finder must not choke when trying to find the roots of
polynomials of the form

${p_{d,\mathbf{e}}}\\!\left(x\right):=\pm\prod_{i=-m}^{m}\left({x-i}\right)^{e_{i}}\hskip
11.38092pt\mathrm{for}\hskip 11.38092ptx\in\mathbb{x}=[m-\underline{\delta},\
m+\overline{\delta}]$ (21)

where $m$ and the exponents $e_{i}$ are integers and

$\underline{\delta},\overline{\delta}\in\left\\{{0,1}\right\\},\hskip
11.38092pt\sum_{i=-m}^{m}e_{i}=d>0,\hskip 11.38092pt\hskip 11.38092ptm>0\hskip
11.38092pt\mathrm{and}\hskip 11.38092pte_{i}\geq 0\hskip
5.69046pt\mathrm{for}\hskip 5.69046pti=-m,\dots,m.$

If $d$ and $m$ are not large then the coefficients of $p_{d,\mathbf{e}}$ can
be computed exactly in floating point arithmetic. For instance, when $d=20$
and $m=5$ the coefficients of $p_{d,\mathbf{e}}$ can be computed exactly with
the usual IEEEE 754 double precision arithmetic. There are

$n_{d,m}:=8\ \binom{2m+d}{d}$ (22)

elements in the family of polynomials in Equation
$\left({\ref{cases}}\right)$. For $m=5$ we have that

$n_{5}=\sum_{d=1}^{20}n_{d,m}\approx 700\ \mathrm{million},$

and for tolerances $\tau_{x}=\tau_{w}=10^{-6}$ and $\tau_{c}=0.001$ with our
code we can test all these $700$ million cases in a couple of days on our
desktop machine, which has an AMD® Ryzen 7 2700x eight-core processor with
64MB of RAM (see the discussion about tolerances in Section 3.)

The family of functions in Equation $\left({\ref{cases}}\right)$ has the
following features, which make it a good test suit:

  1. 1.

By taking $\underline{\delta}=0$ or $\overline{\delta}=0$ we can check how our
code handles roots on the border of the interval $\mathbb{x}$. This situation
is a common source of bugs in root finders.

  2. 2.

Usually some $e_{i}$ will be greater than one, and the polynomials
$p_{d,\mathbf{e}}$ tend to have multiple roots. In fact, they may have roots
with multiplicity as high as $d$, or a couple of roots with multiplicity
$d/2$. These problems are quite challenging and some root finders may return
thousands of candidates intervals for containing such multiple roots.

  3. 3.

It is easy to write an algorithm to generate these polynomials by writing them
in the factored form

${p_{d,\mathbf{e}}}\\!\left(x\right)=\pm q_{1}(-x)x^{e_{0}}q_{2}(x)$ (23)

where the $q_{k}$ are polynomials of the form

${q_{k}}\\!\left(x\right)=\prod_{i=1}^{m}\left({x-i}\right)^{e_{i}}$

with degree at most $d$. For $m=5$ and $d\leq 20$ there are

$n_{q}=\sum_{j=1}^{20}\binom{m+j-1}{j}=53129$

polynomials $q_{k}$ and we can keep them in a table in memory, and let several
threads build the polynomials $p_{d,\mathbf{e}}$ from them using Equation
$\left({\ref{decomp}}\right)$.

  4. 4.

When the coefficients $c_{k}$ of the expanded version of $p_{d,\mathbf{e}}$

${p_{d,\mathbf{e}}}\\!\left(x\right)=\sum_{k=0}^{d}c_{k}x^{k}$ (24)

can be computed exactly, we know the exact roots of the polynomial in Equation
$\left({\ref{expand}}\right)$ and we can use this knowledge to evaluate the
robustness of our code with respect to the inaccuracies in the evaluation of
$p_{d,\mathbf{e}}$ and its derivatives using the expanded form. In other
words, we should test our code using Horner’s method or a Taylor model version
of it to evaluate $p_{d,\mathbf{e}}$ using its expanded form
$\left({\ref{expand}}\right)$, because the evaluation of the expanded form is
usually much less accurate than the evaluation of the product form in Equation
$\left({\ref{cases}}\right)$, and the purpose of testing is to stress our
code.

Once we have chosen the family of functions and the intervals on which will
search for roots, we must decide which tolerances to use for testing. Our
advice is to use both large and small tolerances. Using large tolerances, like
$0.1$, $0.01$ and $0.001$, we can detect gross mistakes, like mistyping. We
may miss some of these gross mistakes if we only use tiny tolerances, because
with tiny tolerances the part of the code affected by mistyping may never be
executed during the tests, or may only cause bugs which are not severe enough
to be detected (we say this based on our personal experience with our own
blunders.) Tests with small tolerances are necessary to access the algorithm’s
accuracy and to check how much time it takes to find accurate roots.

## References

  * [1] Alefeld, G., On the convergence of some interval arithmetic modiﬁcations of Newton’s method, SIAM J. Numer. Anal. 21, 363–372 (1984)
  * [2] Hansen, E., On Solving Systems of Equations Using Interval Arithmetic. Math. Comp. 22, 374–384 (1968)
  * [3] Hansen, E., Interval Forms of Newtons Method, Computing 20, 153–163 (1978)
  * [4] Hansen, E., Walster, G., Global Optimization and Interval Analysis, Second Edition, Revised and Expanded, Marcel Dekker (2004)
  * [5] Kearfott, R. B., and Du, K., The cluster problem in global optimization, the univariate case. Computing Supplement 9, 117–127 (1992)
  * [6] Kearfott, R. B., and Du, K., The cluster problem in multivariate global optimization. Journal of Global Optimization, 5:253–265 (1994)
  * [7] Krawczyk, R. and Neumaier, A., Interval Slopes for Rational Functions and Associated Centered Forms, SIAM J. Numer. Anal. 22, 604–616 (1985)
  * [8] 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)
  * [9] Neumaier, A., Taylor Forms–Use and Limits. Reliable Computing 9, 43–79 (2003)
  * [10] Moore, R. E., Interval Analysis. Prentice-Hall, Englewood Cliffs, NJ. (1966)
  * [11] Nickel, K., On the Newton Method in Interval Analysis. Mathematics Research Center Report 1136, University of Wisconsin (1971)
  * [12] Ratschek, H., Centered forms, SIAM J. Numer. Anal. vol 15, no. 5, 656-662 (1980)
  * [13] 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)