\documentstyle[11pt]{article}
\setlength{\oddsidemargin}{0in}
\setlength{\evensidemargin}{0in}
\setlength{\textwidth}{6.5in}
\setlength{\topmargin}{0in}
\setlength{\textheight}{8.5in}
\newcommand{\bigsqcap}{
\mbox{\Large $ \sqcap$}}
\hbadness=10000

\newtheorem{definition_int}{Definition}
\newenvironment{definition}{\begin{definition_int}\begin{rm}}%
{\end{rm}\end{definition_int}}

\newtheorem{specification_int}{Specification}
\newenvironment{specification}{\begin{specification_int}\begin{rm}}%
{\end{rm}\end{specification_int}}

\newtheorem{implementation_int}{Implementation}
\newenvironment{implementation}{\begin{implementation_int}\begin{rm}}%
{\end{rm}\end{implementation_int}}

\newtheorem{proposition_int}{Proposition}
\newenvironment{proposition}{\begin{proposition_int}\begin{rm}}%
{\end{rm}\end{proposition_int}}

\newtheorem{example_int}{Example}
\newenvironment{example}{\begin{example_int}\begin{rm}}%
{\end{rm}\end{example_int}}

\newtheorem{lemma_int}{Lemma}
\newenvironment{lemma}{\begin{lemma_int}\begin{rm}}%
{\end{rm}\end{lemma_int}}

\newtheorem{remark_int}{Remark}
\newenvironment{remark}{\begin{remark_int}\begin{rm}}%
{\end{rm}\end{remark_int}}

\newenvironment{proof}{\noindent {\bf Proof }}%
{$\Box$
\vspace{.4cm}}

\newtheorem{theorem_int}{Theorem}
\newenvironment{theorem}{\begin{theorem_int}\begin{rm}}%
{\end{rm}\end{theorem_int}}

\newcommand{\separe}[1]{\noindent
\begin{minipage}{1cm}
\setlength{\unitlength}{1cm}
\begin{picture}(.1,.1)(0,0)
        \put(.4,0){\line(1,0){15}}
\end{picture}
\end{minipage}
}

\newcommand{\R}{\Re}
\newcommand{\Rn}{\Re^{n}}
\newcommand{\Ri}{\Re^{\infty}}
\newcommand{\Rp}{\Re^{+}}
\newcommand{\F}{{\cal F}}
\newcommand{\I}{{\cal I}}
\newcommand{\In}{{\cal I}^{n}}
\newcommand{\ap}[1]{\overline{{#1}}}
\newcommand{\floor}[1]{\lfloor {#1} \rfloor}
\newcommand{\ceil}[1]{\lceil {#1} \rceil}
\newcommand{\LEFT}{{\it left}}
\newcommand{\RIGHT}{{\it right}}
\newcommand{\CENTER}{{\it center}}
\newcommand{\fto}{\rightarrow}

\newcommand{\lbra}{{\rm [\hspace{-.02in}[}}
\newcommand{\rbra}{{\rm ]\hspace{-.02in}]}}
\newcommand{\sem}[1]{\lbra {#1} \rbra}
\newcommand{\Snseq}{D_1,\ldots,D _n}
\newcommand{\xnseq}{x_1,\ldots,x_n}
\newcommand{\rnseq}{r_1,\ldots,r_n}
\newcommand{\Xnseq}{x_1,\ldots,x_n}
\newcommand{\Inseq}{I_1,\ldots,I_n}
\newcommand{\IInseq}{I_1,\ldots,I_{i-1},\ap{r_i},I_{i+1},\ldots,I_n}
\newcommand{\Sonseq}{D^0_1,\ldots,D^0_n}
\newcommand{\Ionseq}{I^0_1,\ldots,I^0_n}
\newcommand{\Pmseq}{P_1,\ldots,P_m}
\newcommand{\IE}{{\it IE}}
\newcommand{\NE}{{\it NE}}
\newcommand{\EF}{{\it EF}}
\newcommand{\EE}{{\it EE}}
\newcommand{\TE}{{\it TE}}
\newcommand{\DE}{{\it DE}}
\newcommand{\newton}{{\tt Newton}}
\newcommand{\Bool}{{\it Bool}}
\newcommand{\bo}[1]{{\it box}\{{#1}\}}
\newcommand{\Constraint}{{\it Constraint}}
\newcommand{\Function}{{\it Function}}
\newcommand{\CONSTRAINT}{{\it CONSTRAINT}}
\newcommand{\FUNCTION}{{\it FUNCTION}}
\newcommand{\UDIV}{{\it UDIV}}
\newcommand{\CS}{{\cal S}}
\newcommand{\CON}{{\it cond}}

%____________________________________________________________________
\begin{document}
\Large
\noindent
{\bf Solving Polynomial Systems Using a Branch and Prune Approach}
\normalsize
\\
\\
\parbox{2.1in}{
\begin{center}
{\bf P. Van Hentenryck} \\
Brown University \\
Box 1910 \\
Providence, RI 02912 \\
{\tt pvh@cs.brown.edu} \\
\end{center}
} \ \
\parbox{2.1in}{
\begin{center}
{\bf D. McAllester} \\
MIT AI Lab \\
Technology Square, 545 \\
Cambridge, USA \\
{\tt dam@ai.mit.edu} \\
\end{center}
} \ \
\parbox{2.3in}{
\begin{center}
{\bf D. Kapur} \\
SUNY at Albany \\
Dep. of Computer Science \\
Albany, NY-12222 \\ 
{\tt kapur@cs.albany.edu} 
\end{center}
} 

\section*{Abstract}

This paper presents {\tt Newton}, a branch \& prune algorithm to find
all isolated solutions of a system of polynomial constraints. {\tt
  Newton} can be characterized as a global search method which uses
intervals for numerical correctness and for pruning the search space
early. The pruning in {\tt Newton} consists in enforcing at each node
of the search tree a unique local consistency condition, called
box-consistency, which approximates the notion of arc-consistency
well-known in artificial intelligence. Box-consistency is parametrized
by an interval extension of the constraint and can be instantiated to
produce Hansen-Segupta narrowing operator (used in interval Newton
methods) as well as new operators which are more effective when the
computation is far from a solution.  {\tt Newton} has been evaluated
on a variety of benchmarks from kinematics, chemistry, combustion,
economics, and mechanics. On these benchmarks, it outperforms the
interval Newton methods we are aware of and compares well with
state-of-the-art continuation methods.  Limitations of {\tt Newton}
(e.g. a sensibility to the size of the initial intervals on some
problems) are also discussed. Of particular interest is the
mathematical and programming simplicity of the method.

\section{Introduction}

Many applications in science and engineering (e.g. chemistry,
robotics, economics, mechanics) require finding all isolated solutions
to a system of polynomial constraints over real numbers.  This problem
is difficult due to its inherent computational complexity (i.e., it is
NP-hard) and due to the numerical issues involved to guarantee
correctness (i.e., finding all solutions) and to ensure termination.
Several interesting methods have been proposed in the past for this
task, including two fundamentally different methods: interval Newton
methods (e.g.  \cite{Moore66,Hansen81,Hansen83,Krawczyk69}) and
continuation methods (e.g. \cite{Morgan87,Cools94}). Continuation
methods have been shown to be effective for problems for which the
total degree is not too high, since the number of paths explored
depends on the estimation of the number of solutions. Interval Newton
methods are generally robust but tend to be slow.

The purpose of this paper is to propose and to study a novel algorithm
called {\tt Newton}. From a user standpoint, {\tt Newton} receives as
input a system of polynomial constraints over, say, variables
$x_1,\ldots,x_n$ and a box, i.e., an interval tuple $\langle
I_1,\ldots,I_n \rangle$ specifying the initial range of these
variables; it returns a set of boxes of specified accuracy containing
all solutions.

Operationally, {\tt Newton} is a branch \& prune algorithm which was
inspired by the traditional branch and bound approach used to solve
combinatorial optimization problems. {\tt Newton} uses intervals to
address the two fundamental problems listed above.  Numerical
reliability is obtained by evaluating functions over intervals using
outward rounding (as in interval Newton methods).  The complexity
issue is addressed by using constraints to reduce the intervals early
in the search. The pruning in {\tt Newton} is achieved by enforcing a
unique local consistency condition, called box-consistency, at each
node of the search tree. Box-consistency is an approximation of
arc-consistency, a notion well-known in artificial intelligence
\cite{Mackworth77,Montanari74} and used to solve discrete
combinatorial problems in several systems (e.g.,
\cite{ANNALSOR,book}).  Box-consistency is parametrized by an interval
extension operator for the constraint and can be instantiated to
produce various narrowing operators. In particular, box-consistency on
the Taylor extension of the constraint produces a generalization of
Hansen-Segupta operator \cite{Hansen81} well-known in interval Newton
methods. In addition, box-consistency on the natural extension
produces narrowing operators which are more effective when the
algorithm is not near a solution. {\tt Newton} has the following
properties:
\begin{itemize}

\item {\bf Correctness:} {\tt Newton} finds all isolated solutions to
  the system in the following sense: if $\langle v_1,\ldots,v_n
  \rangle$ is a solution, then {\tt Newton} returns at least one box
  $\langle I_1',\ldots,I_n' \rangle$ such that $v_i \in I_i \;\; (1
  \leq i \leq n)$.  In addition, {\tt Newton} may guarantee the
  existence of a unique solution in some or all the boxes in the
  result.

\item {\bf Termination:} {\tt Newton} always terminates in finite
  time. 

\item {\bf Effectiveness:} {\tt Newton} has been evaluated on a
  variety of benchmarks from kinematics, chemistry, combustion,
  economics, and mechanics. It outperforms the interval Newton methods
  we are aware of and compares well with state-of-the-art continuation
  methods on many problems. Interestingly, {\tt Newton} solves the
  Broyden banded function problem \cite{Hansen81} and Mor\'{e}-Cosnard
  discretization of a nonlinear integral equation \cite{More79} for
  several hundred variables. 

\item {\bf Simplicity and Uniformity:} {\tt Newton} is based on simple
  mathematical results and is easy to use and to implement. It is also
  based on a single concept: box-consistency.
\end{itemize}

The rest of this paper is organized as follows. Section \ref{overview}
gives an overview of the approach. Section \ref{background} contains
the preliminaries. Section \ref{pruning} presents an abstract version
of the branch \& prune algorithm. Section \ref{implementation}
discusses the implementation of the algorithm. Section
\ref{experimental} describes the experimental results. Section
\ref{related} discusses related work and the development of the ideas
presented here. Section \ref{conclusion} concludes the paper.

\section{Overview of The Approach}
\label{overview}

As mentioned, {\tt Newton} is a global search algorithm which solves a
problem by dividing it into subproblems which are solved recursively.
In addition, {\tt Newton} is a branch \& prune algorithm which means that
it is best viewed as an iteration of two steps
\begin{enumerate}
\item pruning the search space;
\item making a nondeterministic choice to generate two subproblems
\end{enumerate}
until one or all solutions to a given problem are found.

The pruning step is responsible to make sure that some local
consistency conditions are satisfied. It consists of reducing the
intervals associated with the variables so that every constraint
appears to be locally consistent. The local consistency condition of
{\tt Newton} is called box-consistency, an approximation of
arc-consistency, a notion well-known in artificial intelligence
\cite{Mackworth77,Montanari74} and used in many systems (e.g.
\cite{book,CCFD,SCREAMER}) to solve discrete combinatorial search
problems.  Informally speaking, a constraint is arc-consistent if for
each value in the range of a variable there exist values in the ranges
of the other variables such that the constraint is satisfied. {\tt
  Newton} approximates arc-consistency which cannot be computed on
real numbers in general.

The pruning step either fails, showing the absence of solution in the
intervals, or succeeds in enforcing the local consistency condition.
Sometimes, local consistency also implies global consistency as in the
case of the Broyden banded function, i.e.
\[
f_i(x_1,\ldots,x_n) = x_i (2 + 5 x_i^2) + 1 - \sum_{j \in J_i} x_j (1 + x_j) \;\;\; (1 \leq i \leq n)
\]
where $J_i = \{ j \ \mid \ j \neq i \ \& \ {\it max}(1,i-5) \leq j
\leq {\it min}(n,i+1)\}$. The pruning step of {\tt Newton} solves this
problem in essentially linear time for initial intervals of the form
$[-10^8,10^8]$ and always proves the existence of a solution in the
final box. However, in general, local consistency does not imply
global consistency either because they are multiple solutions or
simply because the local consistency condition is too weak. Consider
the intersection of a circle and of a parabola:
\[
\left\{ \begin{array}{c}
x_1^2 + x_2^2 = 1 \\
x_1^2 + x_2   = 0
\end{array} \right.
\]
with initial intervals in $[-10^8,10^8]$. The pruning step returns the
intervals
\begin{quote}
$[-1.0000000000012430057,+1.0000000000012430057] $ \\
$[-0.0000000000000000000,+1.0000000000012430057] $
\end{quote}
with 8 digits of accuracy. Branching on $x_1$ directly leads to the
two solutions 
\begin{quote}
$[-0.7861513777574234974,-0.7861513777574231642] $ \\
$[+0.6180339887498946804,+0.6180339887498950136] $ 
\end{quote}
and 
\begin{quote}
$[+0.7861513777574231642,+0.7861513777574233864] $ \\
$[+0.6180339887498946804,+0.6180339887498950136]$.
\end{quote}
Note that, in this case, {\tt Newton} makes the smallest number of
choices to isolate the solutions. To conclude this motivating section,
let us illustrate {\tt Newton} on a larger example which describes the
inverse kinematics of an elbow manipulator \cite{Hong95}:
\[
\left\{ \begin{array}{l}
s_2  c_5  s_6 - s_3  c_5  s_6 - s_4  c_5  s_6 + c_2  c_6 + c_3  c_6 + c_4  c_6 = 0.4077 \\
c_1  c_2  s_5 + c_1  c_3  s_5 + c_1  c_4  s_5 + s_1  c_5 = 1.9115\\
s_2  s_5 + s_3  s_5 + s_4  s_5 = 1.9791\\
c_1  c_2 + c_1  c_3 + c_1  c_4 + c_1  c_2 + c_1  c_3 + c_1  c_2 = 4.0616\\
s_1  c_2 + s_1  c_3 + s_1  c_4 + s_1  c_2 + s_1  c_3 + s_1  c_2 = 1.7172\\
s_2 + s_3 + s_4 + s_2 + s_3 + s_2 = 3.9701\\
s_i^2 + c_i^2 = 1  \;\;\;\;\;\; (1 \leq i \leq 6).
\end{array} \right.
\]
and assumes that the initial intervals are in $[-10^8,10^8]$ again.
The pruning step returns the intervals
\begin{quote}
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[+0.3233666666666665800,+1.0000000000000000000]$  \\
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[-0.0149500000000000189,+1.0000000000000000000]$ \\
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[-0.0209000000000001407,+1.0000000000000000000]$ \\
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[+0.6596999999999998420,+1.0000000000000000000]$ \\
$[-0.7515290480087772896,+0.7515290480087772896]$ \\
$[-1.0000000000000000000,+1.0000000000000000000]$ \\
$[-1.0000000000000000000,+1.0000000000000000000]$ 
\end{quote}
showing already some interesting pruning. After exactly 12 branchings
and in less than a second, {\tt Newton} produces the first box
with a proof of existence of a solution in the box.

\section{Preliminaries}
\label{background}

In this section, we review some basic concepts needed for this paper,
including interval arithmetic and the representation of constraints.
More information on interval arithmetic can be found in many places
(e.g., \cite{Alefeld83,Hansen81,Hansen83,Moore66,Moore79}). Our
definitions are slightly non-standard.

\subsection{Interval Arithmetic}

We consider $\Ri = \R \cup \{-\infty,\infty\}$ the set of real numbers
extended with the two infinity symbols and the natural extension of
the relation $<$ to this set. We also consider a finite subset $\F$ of
$\Ri$ containing $-\infty,\infty,0$. In practice, $\F$ corresponds to
the floating-point numbers used in the implementation. 

\begin{definition}[Interval] An interval $[a,b]$ with $a,b \in \F$ is
the set of real numbers $$\{ r \in \R \ \mid \ a \leq r \leq b \}.$$
The set of intervals is denoded by $\I$ and is ordered by set
inclusion.
\end{definition}

\begin{definition}[Approximation] Let $S$ be a subset of $\R$. The
  approximation of $S$, denoted by $\ap{S}$ or ${\it box\{S\}}$, is the
  smallest interval $I$ such that $S \subseteq I$. We often write $r$
  instead of $\ap{\{r\}}$. We also write $I_1 \uplus I_2$ to denote
  ${\it box}\{I_1 \cup I_2\}$.
\end{definition}

\noindent In the following, we denote real numbers by the letters
$r,v$, $\F$-numbers by the letters $a,b,l,m,u$, intervals by the
letter $I$, real functions by the letters $f,g$ and interval functions
by the letters $F,G$, all possibly subscripted. We use $a^+$ (resp.
$a^-$) to denote the smallest (resp. largest) $\F$-number strictly
greater (resp. smaller) than the $\F$-number $a$. To capture outward
rounding, we use $\ceil{r}$ (resp.  $\floor{r}$) to return the
smallest (resp. largest) $\F$-number greater (resp. smaller) or equal
to the real number $r$. We also use $\vec{I}$ to denote a box $\langle
I_1, \ldots, I_n \rangle$ and $\vec{r}$ to denote a tuple $\langle
r_1, \ldots, r_n \rangle$.  Finally, we use the following notations.
\begin{tabbing}
12345\= \kill
\> $ \LEFT([l,u]) = l$ \\
\> $ \RIGHT([l,u]) = u$ \\
\> $\CENTER([l,u]) = \floor{(l+u)/2}$ 
\end{tabbing}

\noindent The fundamental concept of interval arithmetic is the notion of
interval extension. 

\begin{definition}[Interval Extension] $F: \In \fto \I$ is an interval
extension of $f: \Rn \fto \R$ iff
\[
\forall I_1 \ldots \forall I_n: r_1 \in I_1, \ldots, r_n \in I_n \Rightarrow
f(r_1,\ldots,r_n) \in F(I_1,\ldots,I_n).
\]
An interval relation $C:
  \In \fto {\it Bool}$ is an interval extension of a relation $c: \Rn
  \fto {\it Bool}$ iff
\[
\forall I_1 \ldots \forall I_n: r_1 \in I_1, \ldots, r_n \in I_n
\Rightarrow [c(r_1,\ldots,r_n) \Rightarrow C(I_1,\ldots,I_n)].
\]
\end{definition}

\begin{example}
The interval function $\oplus$ defined as
\[
[a_1,b_1] \oplus [a_2,b_2] = [\floor{a_1 + a_2}, \ceil{b_1 + b_2}]
\]
is an interval extension of addition of real numbers.  The
interval relation $\approx$ defined as
\[
I_1 \approx I_2 \; \Leftrightarrow \; (I_1 \cap I_2 \neq \emptyset)
\]
is an interval extension of the equality relation on real numbers.
\end{example}

\noindent In the following, we assume fixed interval extensions for
the basic real operators $+,-,\times$ and exponentiation (for
instance, the interval extension of $+$ is defined by $\oplus$) and
the basic real relations $=,\geq$. In addition, we overload the real
symbols and use them for their interval extensions. Finally, we denote
relations by the letter $c$ possibly subscripted, interval relations
by the letter $C$ possibly subscripted. Note that constraints and
relations are used as synonyms in this paper.

\subsection{Unions of Intervals}
\label{union}

It is sometimes useful to work with unions of intervals instead of
with single intervals to approximate sets of real numbers, although
unions of intervals are less effective in general.  A typical use of
unions of intervals in interval arithmetic is to take the intersection
of an interval $I$ with the result of an operation of the form $I_1 -
I_2/I_3$. When $0 \notin I_3$, it is easy to compute the resulting
interval by using an interval extension of division. Moreover,
computing $I_1 - I_2/I_3$ with unions of intervals makes it possible
to obtain interesting results even when $0 \in I_3$. The computation
returns either a single interval $I_r$ or two intervals $I_n$ and
$I_p$ and the final result is either $I_r \cap I$ or $(I_n \cap I)
\uplus (I_p \cap I)$. Note that the result of an interval division
$[a,b] / [c,d]$ using unions of intervals where $c \leq 0 \leq d$ and
$c < d$ is specified as follows:
\begin{tabbing}
12345\=1234567890123456789012345667890\=\kill
\> $[\floor{b/c},\infty]$ \> if $b \leq 0$ and $d = 0$ \\
\> $[-\infty,\ceil{b/d}] \cup [\floor{b/c},\infty]$ \> if $b \leq 0$ and
$c < 0 < d$ \\
\> $[-\infty,\ceil{b/d}]$ \> if $b \leq 0$ and $c = 0$ \\
\> $[-\infty,\infty]$ \> if $a < 0 < b$ \\
\> $[-\infty,\ceil{a/c}]$ \> if $a \geq 0$ and $d = 0$ \\
\> $[-\infty,\ceil{a/c}] \cup [\floor{b/c},\infty]$ \> if $a \geq 0$ and
$c < 0 < d$ \\
\> $[\floor{a/d},\infty]$ \> if $a \geq 0$ and $c = 0$.
\end{tabbing}
while $[a,b]/[0,0] = [-\infty,\infty]$.


\subsection{Constraint Representations}

It is well-known that different computer representations of a real
function produce different results when evaluated with floating-point
numbers on a computer. As a consequence, the way constraints are
written may have an impact on the behaviour on the algorithm. For this
reason, a constraint or a function in this paper is considered to be
an expression written in terms of the real functional and relational
symbols. This can be formalized by assuming that all constraints used
in this paper are written in a language whose abstract syntax is
specified by the following grammar: \\ \\ 
%
%
\begin{tabular}{lll}
  $c$&$ \in$&$ \Constraint$\\ $f$&$ \in$&$ \Function$\\ $x_i$&$ \in$&$
  {\it RealVariables} = \{\xnseq\}$\\ $q$&$ \in$&$ {\cal Q}$\\ $n$&$
  \in$&$ {\cal N}$\\ 
\end{tabular} \\ \\

\noindent\begin{tabular}{lll}
$c$&$ ::=$&$ f = 0 \ | \ f \geq 0 $ \\
$f$&$ ::=$&$ q \ | \ x_i \ | \ f + f \ | \ f - f \ | \ f \times f \ |
\ f^n \ | \ (f) $
\end{tabular} \\ 

\noindent It is easy to extend the language to accommodate more
functions (e.g. {\it sin, cos, $\ldots$}). Note that we assume for
simplicity that all constraints are written using a finite (but
arbitrary large) set of variables $\{\xnseq\}$.  Let $e$ be a tuple
$(\rnseq)$ and $e_{\mid i}$ be $r_i$ $(1 \leq i \leq n)$. The
semantics of the language is given by the following semantic functions
whose signatures are
\begin{tabbing}
12345\=\kill
\> ${\cal S}_c: \Constraint \fto \Rn \fto {\it Bool}$ \\
\> ${\cal S}_f: \Function \fto \Rn \fto \R$ 
\end{tabbing}
and whose semantic equations are specified as follows: \\ \\
\begin{tabular}{lll}
${\cal S}_c \sem{f = 0} $& = &$\lambda e. \ {\cal S}_f \sem{f} e \ = \ 0$ \\
${\cal S}_c \sem{f \geq 0} $& = &$\lambda e. \ {\cal S}_f \sem{f} e \ \geq \ 0$ \\
${\cal S}_f \sem{q} $& = &$\lambda e. \ q$\\
${\cal S}_f \sem{x_i}$& = &$\lambda e. \ e_{\mid i}$\\
${\cal S}_f \sem{f_1 + f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e + {\cal S}_f \sem{f_2} e $\\
${\cal S}_f \sem{f_1 - f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e - {\cal S}_f \sem{f_2} e $\\
${\cal S}_f \sem{f_1 \times f_2} $& = &$\lambda e. \ {\cal S}_f \sem{f_1} e \times {\cal S}_f \sem{f_2} e $\\
${\cal S}_f \sem{f^n} $& = &$\lambda e. \ ({\cal S}_f \sem{f} e)^n$\\
${\cal S}_f \sem{(f)} $& = &$\lambda e. \ {\cal S}_f \sem{f} e$
\end{tabular} \\ 

\noindent In the following, we often abuse notation and use $f$ (resp.
$c$) to denote ${\cal S}_f \sem{f}$ (resp. ${\cal S}_c \sem{c}$ ) and
vice-versa. However, the meaning should be clear from the context.  We
also assume that interval functions and interval relations are defined
using a similar abstract syntax and semantics except that real
variables (i.e. $x_i$) are replaced by interval variables (i.e.
$X_i$) and rational numbers by intervals. We use $\CONSTRAINT$ and
$\FUNCTION$ to denote the set of interval constraints and
functions. For simplicity of exposion, we restrict attention to
equations. It is straightforward to generalize our results to
inequalities (see Section \ref{tricks}).


\section{The Branch \& Prune Algorithm}
\label{pruning}

This section describes the branch \& prune algorithm {\tt Newton}.
Section \ref{box} defines box-consistency, the key concept behing our
algorithm. Section \ref{extensions} shows how box-consistency can be
instantiated to produce various pruning operators achieving various
tradeoffs between accuracy and efficiency. Section \ref{cond} defines
a conditioning operator used in {\tt Newton} to improve the
effectiveness of box-consistency. Section \ref{newton-consistency}
specifies the pruning in {\tt Newton}. Section \ref{algo} describes
the algorithm. Recall that we assume that all constraints are defined
over variables $\xnseq$.

\subsection{Box Consistency}
\label{box}

Box-consistency \cite{ILPS-94-CO} is an approximation of
arc-consistency, a notion well-known in artificial intelligence
\cite{Mackworth77} which states a simple local condition on a
constraint $c$ and the set of possible values for each of its
variables, say $D_1,\ldots,D_n$.  Informally speaking, a constraint
$c$ is arc-consistency if none of the $D_i$ can be reduced by using
projections of $c$.

\begin{definition}[Projection Constraint]
  A projection constraint $\langle c, i \rangle$ is the association of
  a constraint $c$ and of an index $i$ ($1 \leq i \leq n$).
  Projection constraints are denoted by the letter $P$, possibly
  subscripted.
\end{definition}

\begin{definition}[Arc-Consistency] A projection constraint $\langle
  c, i\rangle$ is arc-consistent wrt $\langle D_1, \ldots,D_n \rangle$
  iff $ D_i = D_i \cap \{ r_i \ \mid \ \exists r_1 \in D_1, \ldots,
  \exists r_{i-1} \in D_{i-1}, \ldots, \exists r_{i+1} \in D_{i+1},
  \ldots, \exists r_{n} \in D_n: c(r_1,\ldots,r_n) \}.  $ A constraint
  $c$ is arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$ if each
  of its projections is arc-consistent wrt $\langle D_1,\ldots,D_n
  \rangle$. A system of constraints ${\cal S}$ is arc-consistent wrt
  $\langle D_1,\ldots,D_n \rangle$ if each constraint in ${\cal S}$ is
  arc-consistent wrt $\langle D_1,\ldots,D_n \rangle$.
\end{definition}

\noindent Arc-consistency cannot be computed in general when working
with real numbers and polynomial constraints and simple approximations
to capture machine precision are very expensive to compute. For
instance, a simple approximation of arc-consistency consists in
working with intervals and approximating the set computed by
arc-consistency to return an interval, i.e.
\[
I_i = \bo{I_i \cap \{ \ r_i \ \mid \ \exists r_1 \in I_1, \ldots, \exists
r_{i-1} \in I_{i-1}, \ldots, \exists r_{i+1} \in I_{i+1}, \ldots,
\exists r_{n} \in I_n: c(r_1,\ldots,r_n) \ \}}.
\]
This condition, used in systems like \cite{BNR,CLP-BNR}, is easily
enforced on simple constraints such as
\[
x_1 = x_2 + x_3, \;\; x_1 = x_2 - x_3, \;\; x_1 = x_2 \times x_3
\]
but it is also computationally very expensive for complex constraints
with multiple occurrences of the same variables. Moreover, decomposing
complex constraints into simple constraints entails a substantial loss
in pruning, making this approach unpractical on many applications.

The notion of box-consistency introduced in \cite{ILPS-94-CO} is a
coarser approximation of arc-consistency which provides a much better
trade-off between efficiency and pruning. It consists in replacing the
existential quantification in the above condition by the evaluation of
an interval extension of the constraint on the intervals of the
existential variables.  Since different interval extensions will
produce different results, it is necessary and useful to parametrize
the definition of box-consistency with an operator which
maps a constraint into one of its extensions.

\begin{definition}[Extension Operator] An extension operator is a
  function of signature 
\[
\Constraint  \fto \CONSTRAINT
\]
\end{definition}

\begin{definition}[Box-Consistency] 
  An projection constraint $\langle c, i \rangle$ is box-consistent
  wrt $\langle \Inseq \rangle$ and an extension operator \IE{}
  iff
\[
\IE(c)(I_1,\ldots,I_{i-1},[l,l^+],I_{i+1},\ldots,I_n) \;
\wedge \; \IE(c)(I_1,\ldots,I_{i-1},[u^-,u],I_{i+1},\ldots,I_n).
\]
where $l = \LEFT(I_i)$ and $u = \RIGHT(I_i)$.  A constraint is
box-consistent wrt $\vec{I}$ and \IE{} if each of its projections is
box-consistent wrt $\vec{I}$ and \IE{}.  A system of constraints is
box-consistent wrt $\vec{I}$ and \IE{} iff each constraint in the
system is box-consistent wrt $\vec{I}$ and \IE{}.
\end{definition}

\subsection{Interval Extensions for Box Consistency}
\label{extensions}

Box-consistency strongly depends on the extension operator $\IE$ and
different choices of operators can produce very different (often
incomparable) tradeoffs between pruning and computational complexity.
In this section, we consider the three extension operators used in
\newton{}: natural interval extension, expanded interval extension,
and Taylor interval extension.

\subsubsection{Natural Interval Extension}

The easiest extension operator produces the natural interval
extension of a constraint.  Informally speaking, it consists in
replacing each real number by an interval, each real variable by an
interval variable, each real operation by its fixed interval
extension. Formally, it is captured by the following definition.

\begin{definition}[Natural Interval Extension] The natural interval
  extension of a constraint $c$, denoted by $\NE(c)$ and of signature
  $\NE: \Constraint \fto \CONSTRAINT$ is defined inductively as
  follows:
\\ \\
\begin{tabular}{lll}
$\NE \sem{f = 0} $& = &$\NE \sem{f} = \ap{0}$ \\
$\NE \sem{q} $& = &$\ap{q}$\\
$\NE \sem{x_i} $& = &$X_i$\\
$\NE \sem{f_1 + f_2}  $& = &$\NE \sem{f_1}  + \NE
\sem{f_2}  $\\
$\NE \sem{f_1 - f_2} $& = &$\NE \sem{f_1}  - \NE \sem{f_2}  $\\
$\NE \sem{f_1 \times f_2} $& = &$\NE \sem{f_1}   \times \NE \sem{f_2}  $\\
$\NE \sem{f^n}  $& = &$(\NE \sem{f} )^n$\\
$\NE \sem{(f)}  $& = &$\NE \sem{f} $
\end{tabular} \\ 
\end{definition}

The advantage of this extension is that it preserves the way
constraints are written and hence users of the system can choose
constraint representations particularly appropriate for the problem at
hand. A very nice application where this extension is fundamental is the
Mor\'{e}-Cosnard discretization of a nonlinear integral equation (See
Section \ref{more-sec}). 

\subsubsection{Expanded Interval Extension}

The second interval extension used by \newton{} does not preserve the
way constraints are written but uses an expanded form of the
constraints. The key advantage of this extension is that it allows the
algorithm to enforce box-consistency by applying Newton interval
method on univariate real functions. The real functions are derived
from univariate interval constraints obtained by replacing all but one
variable by their intervals. As a consequence, applying
box-consistency will be particularly efficient, although the pruning
may be weaker than for the natural extension due to the dependency
problem of interval arithmetics.\footnote{Note that it is not always
  necessary to go through the expanded form to obtain the above
  property but \newton{} adopts it for simplicity.} The main use of
this extension is to speed up box-consistency on the natural extension
but it may happen that it gives more precision as well.

\begin{definition}[Expanded Form] A constraint $c$ is in expanded form
  iff it is of the form
\[
p_1 + \ldots + p_n \; = \; 0
\]
where each $p_i$ is of the form $q x_1^{m_1} \ldots x_n^{m_n}$ with $q
\in {\cal Q}$ and $m_i \in {\cal N}$ and any pair of products
$(p_i,p_j)$ ($i \neq j$) is such that $p_i$ and $p_j$ have different
sequences of exponents $m_1,\ldots,m_n$. Given a constraint $c$, we
denote its expanded version by $\EF(c)$.\footnote{The expanded version
can easily be turned into a canonical representation for constraints.}
\end{definition}

\begin{definition}[Expanded Interval Extension] The expanded interval
extension of a constraint $c$, denoted by $\EE(c)$ and of signature
$\EE: \Constraint \fto \CONSTRAINT$, is
\[
\NE(\EF(c)).
\]
\end{definition}


\subsubsection{Taylor Interval Extension Generator}

The last interval extension we introduce is based on the taylor
expansion around a point. Instead of presenting a single extension, we
present a generator of extensions, since it will be important to
generate an extension for a box under consideration.  The generator
has signature
\[
 \In  \fto \Constraint \fto \CONSTRAINT
\]
and it assumes that the constraint which it is applied to is of the
form $f \ = \ 0$ where $f$ denotes a function which has continuous
partial derivatives.  Given these assumptions, the key idea behind the
operator is to apply a Taylor expansion of the function around the
center of the box and to bound the rest of the series using the box.

\begin{definition}[Taylor Interval Extension Generator] Let $c$ be a constraint
  $f = 0$, $f$ be a function with continuous partial derivatives,
  $\vec{I}$ be a box $(\Inseq)$, and $c_i$ be the center of $I_i$. The
  Taylor interval extension generator
  $\TE(\vec{I})(c)(X_1,\ldots,X_n)$ is defined
\[
\NE(f)(\ap{c_1},\ldots,\ap{c_n}) + \sum_{i=1}^{n} \DE(f,i)(I_1,\ldots,I_n) \; (X_i - \ap{c_i}) \;\; = \;\; \ap{0}.
\]
where $\DE(f,i)$ returns an interval extension of $\frac{\partial
  f}{\partial x_i}$. 
\end{definition}

\noindent In the current version of our system, the operator $\DE$ of
signature $(\Function \times {\cal N}) \fto \FUNCTION$ is defined as
$\NE(D(f,i))$ where $D: (\Function \times {\cal N}) \fto \FUNCTION$ is
specified as follows: \\ \\ 
\begin{tabular}{lll}
$D \sem{q,i} $& = &$0$\\
$D \sem{x_i,i} $& = &$1$\\
$D \sem{x_j,i} $& = &$0$ $\;\;$ $(i \neq j)$\\
$D \sem{f_1 + f_2,i}  $& = &$D \sem{f_1,i}  + D
\sem{f_2,i}  $\\
$D \sem{f_1 - f_2,i} $& = &$D \sem{f_1,i}  - D \sem{f_2,i}  $\\
$D \sem{f_1 \times f_2,i} $& = &$f_1 \times D \sem{f_2,i} + f_2 \times D \sem{f_1,i}  $\\
$D \sem{f^n,i}  $& = &$n \times f^{n-1} \times D \sem{f,i}$\\
$D \sem{(f),i}  $& = &$D \sem{f,i} $\\
\end{tabular} \\ 

\subsection{Conditioning}
\label{cond}

It is interesting to note that box-consistency on the Taylor interval
extension is closely related to Hansen-Segupta's operator
\cite{Hansen81}, which is an improvement over Krawczyk's operator
\cite{Krawczyk69}.\footnote{Note that the operators of Krawczyk and
  Hansen-Segupta can be used to prove the existence of a solution. We
  will come back to this point in the implementation section.}
 Hansen and Smith \cite{Hansen67} also argued that these
operators are more effective for a system $\{ f_1 = 0, \ldots, f_n = 0
\}$ wrt a box $\langle \Inseq \rangle$ when the interval matrix
\[
M_{ij} = \DE(f_i,j)(I_1,\ldots,I_n) \; (1 \leq i, j \leq n)
\]
is diagonally dominant, i.e., 
\[
mig(M_{i,j}) \geq \sum_{j=1, j\neq i}^{n} mag(M_{i,j})
\]
where 
\[
{\it mig} ([l,u]) = {\it min}(\mid l \mid, \mid u \mid) \mbox{
  and } {\it mag} ([l,u]) = {\it max}(\mid l \mid, \mid u \mid).
\]
They also suggest a conditioning which consists in multiplying the
linear relaxation by a real matrix which is the inverse of the matrix
obtained by taking the center of $M[i,j]$. The resulting system is
generally solved through Gauss-Seidel iterations, giving
Hansen-Segupta's operator.  {\tt Newton} exploits this idea to improve
the effectiveness of box-consistency on the Taylor interval extension.
The conditioning of {\tt Newton} is abstracted by the following
definition.

\begin{definition}[Conditioning] Let $S = \{ f_1 = 0, \ldots, f_n =
  0\}$. A conditioning of $S$ is a system $S' = \{ f_1' = 0, \ldots,
  f_n' = 0\}$ where
\[
f_i' = \sum_{k=1}^{n} a_{ik} f_k
\]
where $a_{ik} \in {\cal Q}$.
\end{definition}

\noindent In its present implementation, {\tt Newton} uses a
conditioning {\it cond}($\{ f_1 = 0, \ldots, f_n = 0\}$,$\vec{I}$)
which returns a system $\{ f_1' = 0, \ldots, f_n' = 0\}$ such that
\[
f_i' = \sum_{k=1}^{n} A_{ik} f_k
\]
where
\[
\begin{array}{lll}
M_{ij} & = & \DE(f_i,j)(I_1,\ldots,I_n) \; (1 \leq i, j \leq n) \\
B_{ij} & = & \CENTER(M_{ij}) \\
A_{ij} & = & \left\{ \begin{array}{ll}
              B^{-1}_{ij} & \mbox{ if $B$ is not singular } \\
              {\it if} \; i = j \; {\it then} \; 1 \; {\it else} \; 0 & \mbox{
                otherwise.}
                    \end{array} \right.
\end{array}
\]
Note that the computation of the inverse of $B$ is obtained by
standard floating-point algorithms.

\subsection{Pruning in {\tt Newton}} 
\label{newton-consistency}

We are now in position to  specify the pruning of {\tt Newton}.

\begin{definition}[Newton Consistency]
  A system of constraints ${\cal S}$ is Newton-consistent wrt
  $\vec{I}$ if $\CS$ is box-consistent wrt $\vec{I}$ and $\NE$, $\CS$
  is box-consistent wrt $\vec{I}$ and $\EE$, and $\CON(\CS,\vec{I})$
  is box-consistent wrt $\vec{I}$ and $\TE(\vec{I})$.
\end{definition}

\noindent The key idea behind {\tt Newton} is to apply
Newton-consistency at each node of the search tree, i.e., {\tt Newton}
reduces the current of values for the variables in such a way that the
constraint system is Newton-consistent wrt the reduced intervals and
no solution is removed.  This is captured by the following definition.

\begin{definition}[Maximal Newton-Consistency] $\vec{I}$
  is maximally Newton-consistent wrt a system of constraints ${\cal
    S}$ and an initial box $\vec{I_0}$ iff $\vec{I}$ is the largest
  box included in $\vec{I_0}$ such that ${\cal S}$ is
  Newton-consistent wrt $\vec{I}$. In the following, we denote  the
  box $\vec{I}$ by {\tt PRUNE}($\CS,\vec{I_0}$)
\end{definition}

\noindent The following soundness result follows from the definition
of interval extensions and the soundness of arc-consistency.

\begin{proposition}[Soundness of Newton-Consistency] 
  Let $\vec{I}$ be maximally Newton-consistent wrt a system of
  constraints ${\cal S}$ and an initial box $\vec{I_0}$ and let
  $\vec{r} \in \vec{I_0}$. If $\vec{r}$ is a solution of ${\cal S}$, then 
  $\vec{r} \in \vec{I}$.
\end{proposition}

\subsection{The Branch and Prune Algorithm {\tt Newton}}
\label{algo}

We are now ready to present the branch and prune algorithm which is
depicted in Figure \ref{mainalgo}. The algorithm makes use of a a
function {\tt BRANCH} is used to split a box along one dimension to
produce two subboxes.  Note a box is empty as soon as one of its
components is empty. The key idea behind the algorithm is extremely
simple. It applies operation {\tt PRUNE} on the initial box. If the
resulting box is empty, there is no solution. If it is small enough,
the box is returned as result. Otherwise, the algorithm splits the box
into two parts and applies the algorithm recursively.  Our
experimental results assume a simple round-robin heuristics to choose
the variable to be split in operation {\tt BRANCH}.

\begin{figure}
\begin{tt}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=\kill
{\bf procedure}  BranchAndPrune(${\cal S}$: Set of Constraint;
$\vec{I_0}: \In$): Set of $\In$; \\
{\bf begin} \\
\>  $\vec{I}$ := PRUNE(${\cal S}$,$\vec{I_0}$); \\
\>  {\bf if} $\neg$ IsEmpty($\vec{I}$) {\bf then} \\
\>\>    {\bf if} IsSmallEnough($\vec{I}$) {\bf then} \\
\>\>\>       {\bf return} $\{\vec{I}\}$ \\
\>\>    {\bf else} \\
\>\>\>        $\langle \vec{I_1}, \vec{I_2} \rangle$ := BRANCH($\vec{I}$);\\
\>\>\>     {\bf return} BranchAndPrune(${\cal S}$,$\vec{I_1}$) $\cup$
BranchAndPrune(${\cal S}$,$\vec{I_2})$ \\
\>\>        {\bf endif} \\
\>  {\bf else} \\
\>\>    {\bf return} $\emptyset$ \\
\>        {\bf endif} \\
{\bf end}  
\end{tabbing}
\end{tt}
\caption{The Branch and Prune Algorithm}
\label{mainalgo}
\end{figure}

\section{Implementation}
\label{implementation}

The purpose of this section is to show how to implement
Newton-consistency for the various interval extensions.  The key idea
behind the implementation is to define narrowing operators for each
projection. Given these operators, it is easy to enforce
maximal Newton-consistency by applying them until no reduction occurs.
We first describe a simple fixpoint algorithm using narrowing
operators. We then show how box-consistency on the interval extensions
can be implemented. Implementation issues are also discussed.

\subsection{The Fixpoint Algorithm for Newton-Consistency}

\begin{figure}
\begin{tt}
\begin{tabbing}
12345\=12345\=12345\=\kill
{\bf procedure} PRUNE({\bf in} ${\cal S}$: Set of Constraint; {\bf inout} $\vec{I}: \In$) \\
{\bf begin} \\
\> {\bf repeat} \\
\>\> $\vec{I_p}$ := $\vec{I}$; \\
\>\> BOX-PRUNE($\CS$,$\CON(\CS,\vec{I})$,$\vec{I}$);  \\
\> {\bf until} $\vec{I} = \vec{I_p}$; \\
{\bf end} \\
\\
{\bf procedure} BOX-PRUNE(${\cal S},{\cal S'}$: Set of Constraint; {\bf inout}
$\vec{I}: \In$) \\
{\bf begin} \\
\> {\bf repeat} \\
\>\> $\vec{I_p}$ := $\vec{I}$; \\
\>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\EE) \ \mid \ c\in {\cal
  S} \; \& \; 1 \leq i \leq n \}$;  \\
\>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\NE) \ \mid \ c\in {\cal
  S} \; \& \; 1 \leq i \leq n \}$;  \\
\>\> $\vec{I}$ := $\bigcap \{ \mbox{{\tt BOX-NARROW}}((c,i),\vec{I},\TE(\vec{I})) \ \mid \ c\in {\cal
  S'} \; \& \; 1 \leq i \leq n \}$;  \\
\> {\bf until} $\vec{I} = \vec{I_p}$; \\
{\bf end} 
\end{tabbing}
\end{tt}
\caption{A Fixpoint Algorithm for Operation {\tt PRUNE}}
\label{fixpoint}
\end{figure}

\noindent Newton-consistency can be enforced easily by using narrowing
operators deduced from the definition of box-consistency. These
operators are used to reduce the interval of a variable using a
projection constraint.

\begin{definition}[Box-Narrowing]
  Let $\langle c,i \rangle$ be a projection constraint, $\langle
  \Inseq \rangle$ be a box, and $\IE$ be a interval extension
  operator. The narrowing operator {\tt BOX-NARROW} is defined as
\begin{tt}
\begin{tabbing}
BOX-NARROW($\langle c, i\rangle$,$\langle \Inseq \rangle$,\IE) = $\langle
I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle$
\end{tabbing}
\end{tt}
where $I$ is defined as the largest set included in $I_i$ such that
$\langle c,i\rangle$ is box-consistent with respect to $\langle I_1,
\ldots, I_{i-1}, I , I_{i+1}, \ldots, I_{n} \rangle$ and $\IE$.
\end{definition}

Maximal Newton-consistency is essentially obtained by applying the
narrowing operators of each projection until no further reduction
occurs. This fixpoint algorithm is depicted in Figure \ref{fixpoint}.
Termination of the algorithm is guaranteed since the set $\F$ is
finite and thus the intervals can only be reduced finitely often.

\subsection{Box-Consistency on the Natural Interval Extension}

\begin{figure}
\begin{tt}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=\kill
{\bf function}  {\tt LNAR}($F,F': \I \rightarrow \I$, $I: \I$): $\I$; \\
{\bf begin} \\
\>  $I$ := $N^*(F,F',I)$; \\
\>  {\bf if} 0 $\in$ $F([\LEFT(I),\LEFT(I)^+])$ {\bf then} \\
\>\>     {\bf return} $I$ \\
\>  {\bf else if} 0 $\notin$ $F(I)$ {\bf then} \\
\>\>     {\bf return} $\emptyset$ \\
\>  {\bf else if} \\
\>\>     $\langle I_1, I_2 \rangle$ := SPLIT($I$); \\
\>\>     {\bf if} ${\tt LNAR}(F,F',I_1) \neq \emptyset$ {\bf then} \\
\>\>\>       {\bf return} ${\tt LNAR}(F,F',I_1) \cup I_2$ \\
\>\>     {\bf else} \\
\>\>\>       {\bf return} ${\tt LNAR}(F,F',I_2)$ \\
\>  {\bf endif} \\
{\bf end;}  
\end{tabbing}
\end{tt}
\caption{The function {\tt LNAR}}
\label{LNAR}
\end{figure}

Box-consistency can often be reduced to two subproblems which,
informally speaking, consist in shrinking the left (resp. the right)
of an interval $I'$ to the leftmost (resp.  rightmost) zero of a
univariate interval function $F$ in $I'$.  Moreover, for
polynomial systems, we will always have at our disposal a function
$F'$ such that, for all $I \subseteq I'$,
\[
0 \in F(I) \; \Rightarrow \; 0 \in N(F,F',I)
\]
where 
\[
N(F,F',I) = I \ \cap \ \ap{\CENTER(I)} - \frac{F(I)}{F'(I)}.
\]
Note that the division and the top-level substraction are evaluated
with unions of intervals and the intersection returns a single
interval (see Section \ref{union}),

The shrinkings can now be obtained easily using a variation of the
interval Newton method.  Figure \ref{LNAR} depicts a simple procedure
{\tt LNAR} to shrink to the leftmost zero (procedure {\tt RNAR} is
defined similarly). It make uses of $N^*(F,F',I) =
\bigcap_{i=0}^{\infty} I_i$ where
\[
\begin{array}{lll}
I_0     & = & I \\
I_{i+1} & = & N(F,F',I_i) \; \; (0 \leq i) \\
\end{array}
\] 
and of a function {\tt SPLIT} which, given an interval $I$, returns
two intervals $I_1$ and $I_2$ such that $I = I_1 \cup I_2$ and
$\LEFT(I_2) \geq \RIGHT(I_1)$.  We are now in position to define the
narrowing operator for the natural interval extension.

\begin{definition}[Narrowing Operator for the Natural Interval
  Extension] \label{nar-nat} The narrowing operator for the natural
  interval extension is defined as follows:
\begin{tt}
\begin{tabbing}
12345\=123\=\kill
BOX-NARROW($\langle f = 0, i \rangle,\langle \Inseq \rangle, \NE$) =
$ \langle I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle$ \\
{\bf where} \\
\> $I$ \>  = {\tt LNAR}$(F,F',{\tt RNAR}(F,F',I_i))$ \\
\> $F$ \> = $\lambda X. \ \NE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$ \\
\> $F'$ \>  = $\lambda X. \ \DE(f,i)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$
\end{tabbing}
\end{tt}
\end{definition}
It is important to note that box-consistency on the natural extension
(and on the expanded extension as well) can be applied even if the
function is not differentiable. It suffices to omit the application of
operator $N^*$ in the functions {\tt LNAR} and {\tt RNAR}.

\subsection{Box-Consistency on the Expanded Interval Extension}

\noindent We now turn to the narrowing operator for the expanded
interval extension. Box-consistency on the expanded interval extension
can be enforced by using the interval function
\[
F = \lambda X. \ \EE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})
\]
and by applying similar techniques as in the previous section.
However, since the function is in expanded form, it is possible to do
better by using an idea from \cite{Hong95}. The key insight is
to sandwitch $F$ exactly between two real functions $f_l$ and $f_u$
defined as follows:
\[
\begin{array}{lll}
f_l &  =  & \lambda x. {\it left}(F(\ap{x})) \\
f_u &  =  & \lambda x. {\it right}(F(\ap{x})).
\end{array}
\]
Box-consistency can now be enforced by searching the leftmost and
rightmost zeros of these two functions using, say, the interval Newton
method. The key advantage of the expanded interval extension is that
it is easy to define the function $f_l$ and $f_u$ constructively. Let
$F$ be of the form 
\[
F = \lambda X. \ I_1 X^{n_1} + \ldots + I_p X^{n_p}.
\] 
Function $f_l$ is defined as
\[
f_l = \lambda x. \ {\it low}(I_1,x,n_1) + \ldots + {\it low}(I_p,x,n_p)
\]
where \\
\[
{\it low}(I,x,n) = \left\{ \begin{array}{ll}
 {\it left}(I) \ x^n & \mbox{ if } x \geq 0 
\vee j \mbox{ is even } \\
  {\it right}(I) \ x^n & \mbox{ otherwise }\\
\end{array}
\right.
\]
Function $f_u$ is defined as
\[
f_u = \lambda x. \ {\it high}(I_1,x,n_1) + \ldots + {\it high}(I_p,x,n_p)
\]
where \\
\[
{\it high}(I,x,n) = \left\{ \begin{array}{ll}
 {\it right}(I) \ x^n & \mbox{ if } x \geq 0 
\vee j \mbox{ is even } \\
  {\it left}(I) \ x^n & \mbox{ otherwise }\\
\end{array}
\right.
\]
It is easy to see that the two definitions of $f_l$ and $f_u$ are
equivalent. The narrowing operator can now be obtained by using Newton
interval method on the above two functions. Some care must be applied
obviously since $f_l$ and $f_u$ are not differentiable at 0. The
method is more efficient than applying the interval Newton method on
the interval function, since intervals have been replaced by numbers
increasing the precision of the Newton operator $N^*$.

\begin{definition}[Narrowing Operator for the Expanded Interval
  Extension] Let $\langle f = 0, i \rangle$ be a projection
  constraint. Let $F$ be $\lambda X. \ 
  \EE(f)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})$, $f_l$ be
  $\lambda x. {\it left}(F(\ap{x}))$, $f_u$ be $\lambda x. {\it
    right}(F(\ap{x}))$. The narrowing operator for the expanded
  interval extension is defined as follows:
\begin{tt}
\begin{tabbing}
123\=123\=\kill
BOX-NARROW($\langle (f = 0, \EE), i \rangle,\langle \Inseq \rangle$) =
$ \langle I_1,\ldots,I_{i-1},[l_i,u_i],I_{i+1},\ldots,I_n \rangle$ \\
{\bf where} \\
\> $l_i$  \> = $\left\{ \begin{array}{ll}
    l & \mbox{ if } 0 \in F([l,l^+]) \\
   {\it left}({\tt LNAR}(F_l,F_l',[l,0]) \cup {\tt
     LNAR}(F_l,F_l',[0,u]) & \mbox{ if } F([l,l^+]) \geq \ap{0} \\
   {\it left}({\tt LNAR}(F_u,F_u',[l,0]) \cup {\tt
     LNAR}(F_u,F_u',[0,u]) & \mbox{ otherwise } \\
                        \end{array}
    \right.$ \\
\> $u_i$  \> = $\left\{ \begin{array}{ll}
    u & \mbox{ if } 0 \in F([u^-,u]) \\
   {\it right}({\tt RNAR}(F_l,F_l',[l,0]) \cup {\tt
     RNAR}(F_l,F_l',[0,u]) & \mbox{ if } F([l,l^+]) \geq \ap{0} \\
   {\it right}({\tt RNAR}(F_u,F_u',[l,0]) \cup {\tt
     RNAR}(F_u,F_u',[0,u]) & \mbox{ otherwise } \\
                        \end{array}
    \right.$ \\
\> $F_l$  \> = $\NE(f_l)$ \\
\> $F_u$  \> = $\NE(f_u)$ \\
\> $F_l'$  \> = $\DE(f_l)$ \\
\> $F_u'$  \> = $\DE(f_u)$ \\
\> $l$    \> = ${\it left}(I_i)$ \\
\> $u$    \> = ${\it right}(I_i)$
\end{tabbing}
\end{tt}
\end{definition}

\subsection{Taylor Interval Extension}

We conclude by presenting the narrowing operator for the Taylor
interval extension.  Box-consistency on the Taylor interval extension
can be enforced by using the interval function
\[
F = \lambda X. \ \TE(f,\langle \Inseq
\rangle)(I_1,\ldots,I_{i-1},X,I_{i+1},\ldots,I_{n})
\]
and by applying a simple implementation of {\tt LNAR} and {\tt RNAR}
where the Newton operator $N^*$ is omitted. However, it is possible to
do better by noticing that the constraint is of the form 
\begin{small}
\[
\NE(f)(\ap{c_1},\ldots,\ap{c_n}) + \sum_{j=1}^{i-1} \DE(f,j)(\vec{I})
(I_j - \ap{c_j}) + \DE(f,j)(\vec{I}) (X_i - \ap{c_i}) + \sum_{j=i+1}^{n}
\DE(f,j)(\vec{I}) (I_j - \ap{c_j}) = \ap{0}
\]
\end{small}
where $c_i = \CENTER(I_i)$ and contains a single variable which can
be isolating to compute box-consistency directly.\footnote{Strictly
  speaking, this technique can increase the precision of
  box-consistency by one floating point number since we remove an
  arithmetic operation.  It could be formalized by associating
  interval extensions with projections.  Doing this would obscure the
  paper without adding any real insight.}

\begin{definition}[Narrowing Operator for the Taylor Interval
  Extension] The narrowing operator for the Taylor interval extension
  is defined as follows:
\[
\mbox{{\tt BOX-NARROW}}(\langle f = 0, i \rangle,\langle \Inseq \rangle,\TE) =  \langle
I_1,\ldots,I_{i-1},I,I_{i+1},\ldots,I_n \rangle 
\]
where
\[
I = I_i \cap \ap{c_i} - \frac{1}{\DE(f,i)(I_1,\ldots,I_n)} 
[\sum_{j=1, j \neq i}^{n} \DE(f,j)(I_1,\ldots,I_n) \; (I_j - \ap{c_i}) - 
\NE(f)(\ap{c_1},\ldots,\ap{c_n})].
\]
and $c_i = \CENTER(I_i)$.
\end{definition}

\noindent Note that the above expression is evaluated using unions of
intervals (see Section \ref{union}).

\subsection{Implementation Issues}
\label{tricks}

We now review some implementation issues which arise in programming the
algorithm.

\paragraph{Priorities} 
In the fixpoint algorithm, it is important for efficiency reasons to
use a priority queue to ensure that projections over the expanded
interval extension be selected first and projections over the Taylor
interval extension be selected last. This is motivated by the fact
that box-consistency over the expanded interval extension can be
enforced efficiently and produces significant pruning when far from a
solution, while box-consistency over the Taylor interval extension
requires the computation of the Jacobian of the systems of constraints
and a preconditioning of the system and is mainly effective when near
a solution. {\tt Newton} also does not enforce box-consistency on the
natural extension when the constraints are given in expanded form and
does not enforce box-consistency on the expanded version whenever it
is believed that the expansion would lose too much precision (e.g. an
expression raised to some power).

\paragraph{Precision} In practice, it is often sufficient to return
intervals whose widths\footnote{The width of $[l,u]$ is $u-l$.} are
below at certain level instead of returning intervals of the form
$[l,l^+]$. It is easy to modify the {\tt BRANCH} operation to split
only intervals whose widths is above the required accuracy. Our system
allows users to specify the accuracy.

\paragraph{Improvement Factor} Box-consistency can sometimes take much
time to remove small parts of the intervals. In these cases, it is
probably more cost-effective to branch. Once again, it is easy to
modify the algorithm to avoid this problem by making sure that the
narrowing operators do not update the intervals unless some
significant reduction has taken place. Since the notion of significant
reduction may be problem-dependent, our system lets users specify the
improvement factor necessary to update an interval in a projection.

\paragraph{Automatic Differentiation} As mentioned, our algorithm
takes a very simple approach to obtain partial derivatives, i.e., no
effort is spent in factoring common expressions to reduce the
dependency problem of interval arithmetic. The main reason comes from
the fact that we are using automatic differentiation \cite{Rall81} to
evaluate the derivatives together with the functions. This choice may
be reconsidered in a further version of the system.

\paragraph{Existence Proof of a Solution} As we mentioned earlier,
box-consistency on a system of equations can be viewed as a
generalization of Hansen-Segupta operator used in interval Newton
methods. This operator can be used to prove the existence of
solutions. Our algorithm is organized to exploit this property so that
it is able to report this information to users. See
\cite{Krawczyk86} for more discussion on this topic.

\paragraph{Inequalities} It is simple to generalize the above
algorithms for inequalities. In general, it suffices to test if the
inequality is satisfied at the end points of the interval. If it is
not, then the problem reduces once again to finding the leftmost and/or
rightmost zeros. 

\section{Experimental Results}
\label{experimental}

This section reports experimental results of {\tt Newton} on a variety
of standard benchmarks. The benchmarks were taken from papers on
numerical analysis \cite{More79}, interval analysis
\cite{Hansen81,Hong95,Moore77a}, and continuation methods
\cite{Cools94,Morgan87,Morgan-kin,Morgan90}. We also compare {\tt
  Newton} with a traditional interval method using Hansen-Segupta's
operator, range testing, and branching. This method uses the same
implementation technology as {\tt Newton} and is denoted by {\tt HRB}
in the following.\footnote{Some interval methods such as
  \cite{Hansen83} are more sophisticated than {\tt HRB} but the
  sophistication aims at speeding up the computation near a solution.
  Our main contribution is completely orthogonal and aims at speeding
  up the computation when far from a solution and hence comparing it
  to {\tt HRB} is meaningful.} Finally, we compare {\tt Newton} with a
state-of-the-art continuation method \cite{Cools94}, denoted by {\tt
  CONT} in the following. Note that all results given in this section
were obtained by running {\tt Newton} on a {\tt Sun Sparc 10}
workstation to obtain all solutions.  In addition, the final intervals
must have widths smaller than $10^{-8}$ and {\tt Newton} always uses
an improvement factor of 10\%. The results are summarized in Table
\ref{summary}. For each benchmark, we give the number of variables
($n$), the total degree of the system ($d$), the initial range for the
variables, and the results of each method in seconds. Note that the
times for the continuation method are on a {\tt DEC 5000/200}.  A
space in a column means that the result is not available for the
method. A question mark means that the method does not terminate in a
reasonable time ($>$ 1 hour).  The rest of the section describes each
benchmark and the results in much more detail.

\begin{table}
\small
\centering
\begin{tabular}{||l|c|c|c||r|r|r||}
\hline
\hline
Benchmarks &  $v$  & $d$  & range &  {\tt Newton}  & {\tt HRB}  & {\tt CONT} \\
\hline
Broyden    &  10   & $3^{10}$ & [-1,1] &  1.65 & 18.23 &   \\
Broyden    &  20   & $3^{20}$  & [-1,1] & 4.25 & ?     &  \\
Broyden    & 320   & $3^{320}$ & [-1,1] & 113.71 & ?   & \\
Broyden    & 320   & $3^{320}$ & $[-10^8,10^8]$ & 143.40 & ? &  \\
Mor\'{e}-Cosnard & 20   & $3^{20}$ & $[-4,5]$ & 24.49  & 968.25 &  \\
Mor\'{e}-Cosnard & 40   & $3^{40}$ & $[-4,5]$ & 192.81 & ?      & \\
Mor\'{e}-Cosnard & 80   & $3^{80}$ & $[-4,5]$ & 1752.64 & ?   & \\
Mor\'{e}-Cosnard & 80   & $3^{80}$ & $[-10^8,0]$ & 1735.09 & ? & \\
i1         & 10    & $3^{10}$ & [-2,2] & 0.06 & 14.28 & \\
i2         & 20    & $3^{20}$ & [-1,2] & 0.30 & 1821.23 & \\
i3         & 20    & $3^{20}$ & [-2,2] & 0.31 & 5640.80 & \\
i4         & 10    & $6^{10}$ & [-1,1] & 73.94 & 445.28 & \\
i5         & 10    & $11^{10}$ & [-1,1] & 0.08 & 33.58 & \\
kin1       & 12    & 4608     &  $[-10^8,10^8]$ & 14.24 & 1630.08 & \\
kin2       & 8     & 256      &  $[-10^8,10^8]$ & 353.06 & 4730.34 & 35.61 \\
eco        & 4     & 18       &  $[-10^8,10^8]$ & 0.60 &  & 1.13 \\
eco        & 5     & 54       &  $[-10^8,10^8]$ & 3.35 &  & 5.87 \\
eco        & 6     & 162      &  $[-10^8,10^8]$ & 22.53 &  & 50.18 \\
eco        & 7     & 486      &  $[-10^8,10^8]$ & 127.65 &  & 991.45 \\
eco        & 8     & 1458     &  $[-10^8,10^8]$ & 915.24 &  &  \\
eco        & 9     & 4374     &  $[-10^8,10^8]$ & 8600.28 &  &  \\
combustion & 10    & 96       &  $[-10^8,10^8]$ & 9.94   &   & 57.40 \\
chemistry  & 5     & 108      &  $[0,10^8]$ & 6.32   &   & 56.55 \\
neuro      & 6     & 1024     &  $[-10,10]$ & 0.91   &   & 5.02 \\
neuro      & 6     & 1024     &  $[-1000,1000]$ & 172.71 &  & 5.02 \\
\hline
\end{tabular}
\caption{Summary of the Experimental Results}
\label{summary}
\end{table}



\subsection{Broyden Banded Functions}

This is a traditional benchmark of interval techniques and was used
for instance in \cite{Hansen83}. It consists in finding the
zeros of the functions
\[
{\textstyle f_i(x_1,\ldots,x_n) = x_i (2 + 5 x_i^2) + 1 - \sum_{j \in J_i} x_j (1 + x_j) \;\;\; (1 \leq i \leq n)}
\]
where $J_i = \{ j \ \mid \ j \neq i \ \& \ {\it max}(1,i-5) \leq j
\leq {\it min}(n,i+1)\}$. One of the interesting features of this
benchmark is that it is easy to scale up to an arbitrary dimension and
hence provides a good basis to compare various methods.  Table
\ref{broyden-newton} reports the results of our algorithm for various
sizes assuming initial intervals $[-1,1]$. For each size, we report
the CPU times in seconds, the growth of the CPU time, the number of
branch operations {\it branching}, the number of narrowings on the
various extensions {\it na-ne, na-ee, na-te}, the total number of
narrowings {\it na-tot}, the number of function evaluations (including
evaluation of derivatives) for each of the extensions {\it fe-ne,
  fe-ee, fe-te} and the total number of function evaluations {\it
fe-tot}.  We also indicate
%the growth of the function evaluations {\it fe-grow} and 
the number of preconditionings by {\it pr-con} and whether the algorithm
can prove the existence of the solutions in the resulting intervals
by {\it proof}.

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r|r|r||}
\hline
\hline
               &  5   &  10  &  20   &   40  &   80  &    160  &   320 \\
\hline
{\it CPU time}    & 0.20 & 1.65 & 4.25 & 9.79  & 22.13  & 48.30   & 113.71 \\
{\it growth}   &      & 8.25 & 2.57 & 2.30  & 2.26   & 2.18    & 2.35  \\
{\it branching}   & 0    & 0    & 0     & 0     & 0     & 0       & 0 \\
\hline 
{\it na-ne}    & 57   & 260  & 661   & 1607  & 4351  & 8096    & 17126 \\
{\it na-ee}    & 1226 & 8334 & 21236 & 48540 & 102797& 206926  & 414798 \\
{\it na-te}    & 35   & 110  & 260   & 560   & 1200  & 2400    & 4480 \\
{\it na-tot}   & 1318 & 8704 & 22157 & 50707 & 108348& 217422  & 436404 \\
\hline   
{\it fe-ne}    & 81   & 1828 & 2943  & 5103  & 11993 & 20431   & 42125 \\
{\it fe-ee}    & 3462 & 21518& 53722 & 121984& 257398& 517640  & 1036210 \\
{\it fe-te}    & 95   & 320  & 920   & 2720  & 8800  & 30400   & 111360 \\
{\it fe-tot}   & 3638 & 23666& 57585 & 129807& 278191& 568471  & 1189695 \\
\hline 
{\it pr-con}   & 0    & 0    & 0     & 0     & 0     &  0      & 0 \\
\hline 
{\it proof}    & yes  & yes  & yes   & yes   & yes   &  yes    & yes \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Broyden Banded functions with initial
  intervals $[-1,1]$}
\label{broyden-newton}
\end{table}

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r|r|r||}
\hline
\hline
             &  5   &  10  &  20   &   40  &   80  &    160  &   320 \\
\hline
{\it CPU time}  & 0.31 & 2.15 & 5.09  & 12.49 & 27.69 &  61.60  & 143.40 \\
{\it growth} &      & 6.93 & 2.36 &  2.45  & 2.21  &  2.22   &  2.32  \\
{\it branching} & 0    & 0    &  0    &  0    &  0    &  0      &  0 \\
\hline 
{\it pr-con} & 0    & 0    &  0    &  0    &  0    &  0      &  0 \\
\hline 
{\it proof}  & yes  & yes  &  yes  &  yes  &  yes  &  yes    &  yes  \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Broyden Banded functions with initial
  intervals $[-10^8,10^8]$}
\label{broyden-newton-W}
\end{table}



The results indicate that {\tt Newton} solves the problem using only
constraint propagation: no branching is needed. In addition, the
growth of the computation times is very low and indicates that {\tt
  Newton} is essentially linear and can thus solve very large
instances of this problem. Finally, {\tt Newton} proves the existence
of a solution in the final intervals.  To our knowledge, no other
algorithm has all these functionalities.  Table \ref{broyden-newton-W}
shows the same results when the initial intervals are $[-10^8,10^8]$.
They indicate that the CPU time increases only slightly in this
problem when the initial intervals become substantially larger.  It is
interesting to note that substantial pruning is obtained by
box-consistency on the natural and expanded extensions alone. For
$n=10$, maximal box-consistency on these two extensions produces the
intervals
\begin{quote}
$[-0.4283028737061274627,-0.4283028534683728794]$ \\
$[-0.4765964317901201786,-0.4765964169224605195]$ \\
$[-0.5196524683730758821,-0.5196524589206473754]$ \\
$[-0.5580993358758108425,-0.5580993137885511545]$ \\
$[-0.5925061654931400579,-0.5925061481657747375]$ \\
$[-0.6245036923913307448,-0.6245036720076052594]$ \\
$[-0.6232394806883442274,-0.6232394621928379896]$ \\
$[-0.6213938520278742273,-0.6213938315652728361]$ \\
$[-0.6204536054436834425,-0.6204535878744913413]$ \\
$[-0.5864692773020701023,-0.5864692641387999616]$ 
\end{quote}
which have widths lower than $10^{-6}$. Note that the Hansen-Segupta's
operator alone does not produce any pruning initially and returns the
initial intervals whether they be of the form $[-10^8,+10^8]$ or
$[-1,1]$. This indicates that box-consistency on the natural and
expanded interval extensions are particularly effective when far from
a solution while box-consistency on the Taylor extension (and the
Hansen-Segupta's operator) is effective when near a solution.

It is also interesting to stress the importance of box-consistency on
the natural extension in this example to reduce the growth factor.
Without it, the algorithm takes about 48 and 440 seconds instead of 27
and 61 for {\tt Newton} for $n = 80$ and $n = 160$, since the expanded
interval extension loses precision due to the dependency problem.

Finally, it is interesting to compare {\tt Newton} with traditional
interval Newton methods. {\tt HRB} takes 0.34 seconds on $n=5$ with 18
branchings, about 18 seconds for $n=10$ with about 300 branchings, and
does not return after more than an hour on $n=20$.

\subsection{Discretization of a Nonlinear Integral Equation} 
\label{more-sec}

This example comes from \cite{More79} and is also a standard benchmark
for nonlinear equation solving . It consists in finding the root of
the functions ${\textstyle f_k(x_1,\ldots,x_m)}$ $( 1 \leq k \leq m)$
defined as
\[
x_k + \frac{1}{2(m+1)} [(1 - t_k) \sum_{j=1}^{k} t_j (x_j + t_j + 1)^3 + t_k \sum_{j=k+1}^{m} (1 -t_j)(x_j + t_j + 1)^3 ]
\]
where $t_j = jh$ and $h = 1/(m+1)$. These functions come from the
discretization of a nonlinear integral equation, giving a constraint
system denser than the sparse constraint system for the Broyden banded
functions.  The variables $x_i$ were given initial domains $[-4,5]$ as
in \cite{Rokne88} and the computation results are given in Table
\ref{more-newton}.
\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
               &  5   &  10  &  20   &   40   &   80   \\
\hline
{\it CPU time} & 0.75 & 4.07 & 24.49 & 192.81 & 1752.64 \\
{\it growth}   &      & 5.42 & 6.02  &   7.87 & 9.08    \\
{\it branching}   & 0    & 0    & 0     & 0      & 0       \\
\hline 
{\it na-ne}    & 3663 & 12616& 46555 & 213949 & 1236532 \\
{\it na-te}    & 104  & 255  & 505   & 1005   & 2807    \\
{\it na-tot}   & 3767 & 12871& 47060 & 214954 & 1239339 \\
\hline   
{\it fe-ne}    & 8775 & 31837& 107977& 466595 & 2586907 \\
{\it fe-te}    &  884 & 3111 & 11211 & 42411  & 166415  \\
{\it fe-tot}   & 9659 & 34948& 119188& 509006 & 2753322 \\
{\it fe-grow}  &      & 3.62 & 3.41  & 4.27   & 5.40    \\
\hline 
{\it pr-con}   & 1    & 1    & 1     & 1      & 1   \\
\hline 
{\it proof}    & yes  & no   & no    & no     & no \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Mor\'{e}-Cosnard nonlinear integral
  Equation with initial intervals in $[-4,5]$}
\label{more-newton}
\end{table}

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
                &  5   &  10  &  20   &   40  &   80  \\
\hline
{\it CPU time}  & 0.70 & 3.82 & 20.81 & 189.94 &  1735.09 \\
{\it growth}    &      & 5.45 & 5.44  &  9.12  &  9.13 \\
{\it branching}    & 0    & 0    &  0    &  0    &  0     \\
\hline 
{\it pr-con}    & 0    & 0    &  0    &  0    &  0      \\
\hline 
{\it proof}     & yes  & no   &  no   &  no   &  no     \\
\hline
\hline 
\end{tabular}
\caption{{\tt Newton} on the Mor\'{e}-Cosnard nonlinear integral
  Equation with initial intervals in $[-10^8,0]$}
\label{more-newton-W}
\end{table}

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
               &  5   &  10  &  20   &   40   &   80   \\
\hline
{\it CPU time} & 0.66 & 7.76 & 968.25     & ?      & ? \\
{\it branching}   &  5   & 24   & 508     & ?      & ?       \\
\hline 
{\it fe-tot}   & 3709 & 20194 & 1285764     & ?      & ? \\
\hline 
{\it pr-con}   & 7    & 32    & 667     & ?      & ?  \\
\hline 
{\it proof}    & yes  & no   & no    & ?      & ? \\
\hline
\hline 
\end{tabular}
\caption{The {\tt HRB} algorithm on the Mor\'{e}-Cosnard
 nonlinear integral Equation with initial intervals in $[-4,5]$}
\label{more-hansen}
\end{table}

Once again, it is interesting to note that {\tt Newton} is completely
deterministic on this problem, i.e. it does not do any branching. {\tt
  Newton} is probably cubic in the number of variables for this
problem.  It is important to point out the critical role of
box-consistency on the natural extension to solve this problem
efficiently. {\tt Newton} without the natural extension would not be
deterministic and would slow down exponentially, since box-consistency
on the expanded extension loses too much precision due to the
dependency problem (multiple occurrences of the same variable) and
box-consistency on the Taylor interval extension is not helpful
initially. Once again, we observe that box-consistency over the
natural extension is helpful when far from a solution while
box-consistency on the Taylor extension is useful to terminate the
search quickly. Table \ref{more-newton-W} gives the result for the
initial intervals of size $[-10^8,0]$, which shows that the algorithm
continues to perform well in this case. Finally, Table
\ref{more-hansen} gives the results for the {\tt HRB} algorithm on
this problem. Once again, {\tt Newton} outperforms the {\tt HRB}
method substantially.

\subsection{Interval Arithmetic Benchmarks}

This section considers standard benchmarks from interval arithmetic
papers \cite{Moore79,Hong95}. Benchmark {\tt i1} is the following set
of equations
\begin{small}
\[
\left\{ \begin{array}{l}
0 = x_{1} -  0.25428722 - 0.18324757 \ x_{4} \ x_{3}  \ x_{9} \\
0 = x_{2} -  0.37842197 - 0.16275449 \ x_{1} \ x_{10}  \ x_{6} \\
0 = x_{3} -  0.27162577 - 0.16955071 \ x_{1} \ x_{2}  \ x_{10} \\
0 = x_{4} -  0.19807914 - 0.15585316 \ x_{7} \ x_{1}  \ x_{6} \\
0 = x_{5} -  0.44166728 - 0.19950920 \ x_{7} \ x_{6}  \ x_{3} \\
0 = x_{6} -  0.14654113 - 0.18922793 \ x_{8} \ x_{5}  \ x_{10} \\
0 = x_{7} -  0.42937161 - 0.21180486 \ x_{2} \ x_{5}  \ x_{8} \\
0 = x_{8} -  0.07056438 - 0.17081208 \ x_{1} \ x_{7}  \ x_{6} \\
0 = x_{9} -  0.34504906 - 0.19612740 \ x_{10} \ x_{6}  \ x_{8} \\
0 = x_{10} - 0.42651102 - 0.21466544 \ x_{4} \ x_{8}  \ x_{1} 
\end{array} \right.
\]
\end{small}
with initial intervals $[-2,2]$. Benchmark {\tt i2} is the 
set of equations
\begin{small}
\[
\left\{ \begin{array}{l}
0 = x_{1}  - 0.24863995 - 0.19594124 \ x_{7} \ x_{10}  \ x_{16} \\
0 = x_{2}  - 0.87528587 - 0.05612619 \ x_{18} \ x_{8}  \ x_{11} \\
0 = x_{3}  - 0.23939835 - 0.20177810 \ x_{10} \ x_{7}  \ x_{11} \\
0 = x_{4}  - 0.47620128 - 0.16497518 \ x_{12} \ x_{15}  \ x_{1} \\
0 = x_{5}  - 0.24711044 - 0.20198178 \ x_{8} \ x_{9}  \ x_{16} \\
0 = x_{6}  - 0.33565227 - 0.15724045 \ x_{16} \ x_{18}  \ x_{11} \\
0 = x_{7}  - 0.13128974 - 0.12384342 \ x_{12} \ x_{13}  \ x_{15} \\
0 = x_{8}  - 0.45937304 - 0.18180253 \ x_{19} \ x_{15}  \ x_{18} \\
0 = x_{9}  - 0.46896600 - 0.21241045 \ x_{13} \ x_{2}  \ x_{17} \\
0 = x_{10} - 0.57596835 - 0.16522613 \ x_{12} \ x_{9}  \ x_{13} \\
0 = x_{11} - 0.56896263 - 0.17221383 \ x_{16} \ x_{17}  \ x_{8} \\
0 = x_{12} - 0.70561396 - 0.23556251 \ x_{14} \ x_{11}  \ x_{4} \\
0 = x_{13} - 0.59642512 - 0.24475135 \ x_{7} \ x_{16}  \ x_{20} \\
0 = x_{14} - 0.46588640 - 0.21790395 \ x_{13} \ x_{3}  \ x_{10} \\
0 = x_{15} - 0.10607114 - 0.20920602 \ x_{1} \ x_{9}  \ x_{10} \\
0 = x_{16} - 0.26516898 - 0.21037773 \ x_{4} \ x_{19}  \ x_{9} \\
0 = x_{17} - 0.20436664 - 0.19838792 \ x_{20} \ x_{10}  \ x_{13} \\
0 = x_{18} - 0.56003141 - 0.18114505 \ x_{6} \ x_{13}  \ x_{8} \\
0 = x_{19} - 0.92894617 - 0.04417537 \ x_{7} \ x_{13}  \ x_{16} \\
0 = x_{20} - 0.57001682 - 0.17949149 \ x_{1} \ x_{3}  \ x_{11}
\end{array} \right.
\]
\end{small}

\noindent with initial intervals $[-1,2]$. Benchmark {\tt i3} has the same set
of equations as {\tt i2} but has initial intervals $[-2,2]$. Benchmark
{\tt i4} has the set of equations
\begin{small}
\[
\left\{ \begin{array}{l}
0 = x_1^2 - 0.25428722 - 0.18324757 \  x_4^2 \  x_3^2 \  x_9^2 \\
0 = x_2^2 - 0.37842197 - 0.16275449 \  x_1^2 \  x_{10}^2 \  x_6^2 \\
0 = x_3^2 - 0.27162577 - 0.16955071 \   x_1^2 \   x_2^2 \   x_{10}^2 \\
0 = x_4^2 - 0.19807914 - 0.15585316 \   x_7^2 \   x_1^2 \   x_6^2 \\
0 = x_5^2 - 0.44166728 - 0.19950920 \   x_7^2 \   x_6^2 \   x_3^2 \\
0 = x_6^2 - 0.14654113 - 0.18922793 \   x_8^2 \   x_5^2 \   x_{10}^2 \\
0 = x_7^2 - 0.42937161 - 0.21180486 \   x_2^2 \   x_5^2 \   x_8^2 \\
0 = x_8^2 - 0.07056438 - 0.17081208 \   x_1^2 \   x_7^2 \   x_6^2 \\
0 = x_9^2 - 0.34504906 - 0.19612740 \   x_{10}^2 \   x_6^2 \   x_8^2 \\
0 = x_{10}^2 - 0.42651102 - 0.21466544 \   x_4^2 \   x_8^2 \   x_1^2
\end{array} \right.
\]
\end{small}

\noindent and initial intervals $[-1,1]$. The number of solutions must be
a multiple of 1024. Benchmark {\tt i5} has the following set of
equations
\begin{small}
\[
\left\{ \begin{array}{l}
0 = x_1 - 0.25428722 - 0.18324757 \  x_4^3 \   x_3^3 \   x_9^3 \ + \
x_3^4 \  
x_9^7 \\
0 = x_2 - 0.37842197 - 0.16275449 \  x_1^3 \   x_{10}^3 \   x_6^3 \ +
\ x_{10}^4 \  
x_6^7 \\
0 = x_3 - 0.27162577 - 0.16955071 \  x_1^3 \   x_2^3 \   x_{10}^3 \ +
\ x_2^4 \  
x_{10}^7 \\
0 = x_4 - 0.19807914 - 0.15585316 \  x_7^3 \   x_1^3 \   x_6^3 \ + \
x_1^4 \  
x_6^7 \\
0 = x_5 - 0.44166728 - 0.19950920 \  x_7^3 \   x_6^3 \   x_3^3 \ + \
x_6^4 \  
x_3^7 \\
0 = x_6 - 0.14654113 - 0.18922793 \  x_8^3 \   x_5^3 \   x_{10}^3 \ +
\ x_5^4 \  
x_{10}^7 \\
0 = x_7 - 0.42937161 - 0.21180486 \  x_2^3 \   x_5^3 \   x_8^3 \ + \
x_5^4 \  
x_8^7 \\
0 = x_8 - 0.07056438 - 0.17081208 \  x_1^3 \   x_7^3 \   x_6^3 \ + \
x_7^4 \  
x_6^7 \\
0 = x_9 - 0.34504906 - 0.19612740 \  x_{10}^3 \   x_6^3 \   x_8^3 \ +
\ x_6^4 \  
x_8^7 \\
0 = x_{10} - 0.42651102 - 0.21466544 \  x_4^3 \   x_8^3 \   x_1^3 \ +
\ x_8^4 \  x_1^7
\end{array} \right.
\]
\end{small}

\noindent and initial intervals $[-1,1]$. 

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
               &{\tt i1} &{\tt i2} &{\tt i3}&{\tt i4}&{\tt i5}   \\
\hline
{\it CPU time} &  0.06   & 0.30    & 0.31   & 73.94   & 0.08 \\
{\it branching}   &  0      & 0       & 0      & 1023    & 0 \\
\hline 
{\it na-ee}    &  625    & 2107    & 1949   & 290776  & 381 \\
{\it na-te}    &  0      & 80      & 80     & 37930   & 21 \\
{\it na-tot}   &  625    & 2187    & 2029   & 328706  & 401 \\
\hline   
{\it fe-ee}    &  1760   & 5698    & 5318   & 752220  & 992 \\
{\it fe-te}    &  0      & 560     & 560    & 269760  & 140 \\
{\it fe-tot}   &  1760   & 6258    & 5878   & 1021980 & 1132 \\
\hline 
{\it pr-con}   &  0      & 1       & 1      & 1939    & 1   \\
\hline 
{\it proof}    &  yes     & yes     & yes    & yes    & yes \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Traditional Interval Arithmetic Benchmarks}
\label{trad-newton}
\end{table}

{\tt Newton} solves all the problems with one solution without
branching and solves the problem having 1024 solutions with 1023
branchings. Note also that box-consistency on the expanded extension
solves benchmark {\tt i1} alone. The results once again confirm our
observation on when the various extensions are useful.  Closely
related results were observed in \cite{Hong95} on these benchmarks
(see the related work section for a more detailed comparison) but our
algorithm is in general about 4 times faster (assuming similar
machines) and does not do any branching on {\tt i5}. Table
\ref{trad-hansen} also describes the results for the traditional
interval arithmetic method.  The importance of box-consistency on the
expanded extension can easily be seen from these results. Note also
that {\tt Newton} (and interval methods) can prove the existence of a
solution in the final intervals for all these problems.

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
               &{\tt i1} &{\tt i2} &{\tt i3}&{\tt i4}&{\tt i5}   \\
\hline
{\it CPU time} &  14.28  & 1821.23 & 5640.80& 445.28 &  33.58     \\
{\it branching}   &  498    & 9031    & 36933  & 11263 &  1173     \\
\hline 
{\it fe-tot}   &  77380  & 6441640 & 19979025 &2554066 &  154948       \\
\hline 
{\it pr-con}   &  586    & 12817   & 42193 & 14335  &  1211    \\
\hline 
{\it proof}    &  yes    & yes     & yes       & yes    &  yes      \\
\hline
\end{tabular}
\caption{{\tt HRB} on the Traditional Interval Arithmetic Benchmarks}
\label{trad-hansen}
\end{table}

\subsection{Kinematics Applications}

We now describe the performance of {\tt Newton} on two kinematics
examples.  Application {\tt kin1} comes from robotics and describes
the inverse kinematics of an elbow manipulator \cite{Hong95}.  It
consists of a sparse system with 12 variables and the set of equations
is as follows:
\[
\left\{ \begin{array}{l}
s_2  c_5  s_6 - s_3  c_5  s_6 - s_4  c_5  s_6 + c_2  c_6 + c_3  c_6 + c_4  c_6 = 0.4077 \\
c_1  c_2  s_5 + c_1  c_3  s_5 + c_1  c_4  s_5 + s_1  c_5 = 1.9115\\
s_2  s_5 + s_3  s_5 + s_4  s_5 = 1.9791\\
c_1  c_2 + c_1  c_3 + c_1  c_4 + c_1  c_2 + c_1  c_3 + c_1  c_2 = 4.0616\\
s_1  c_2 + s_1  c_3 + s_1  c_4 + s_1  c_2 + s_1  c_3 + s_1  c_2 = 1.7172\\
s_2 + s_3 + s_4 + s_2 + s_3 + s_2 = 3.9701\\
s_i^2 + c_i^2 = 1  \;\;\;\;\;\; (1 \leq i \leq 6).
\end{array} \right.
\]

The second benchmark, denoted by {\tt kin2}, is from \cite{Morgan-kin}
and describes the inverse position problem for a six-revolute-joint
problem in mechanics. The equations which describe a denser constraint
system are as follows:
\[
\left\{ \begin{array}{l}
x_i^2 + x_{i+1}^2 - 1 = 0 \;\;\;\;\; (1 \leq i \leq 4) \\
a_{1i} x_1 x_3 + a_{2i} x_1 x_4 + a_{3i} x_2 x_3 + a_{4i} x_2 x_4 
+ a_{5i} x_5 x_7 + a_{6i} x_5 x_8 + a_{7i} x_6 x_7 + a_{8i} x_6 x_8 \\
a_{9i} x_1 + a_{10i} x_2 + a_{11i} x_3 + a_{12i} x_4 + a_{13i} x_5
a_{14i} x_6 + a_{15i} x_7 + a_{16i} x_8 + a_{17i} = 0  \;\; (1 \leq i
\leq 4)
\end{array} \right.
\]
where the coefficients $a_{ki}$ are given in table \ref{kin-coef}. In
both examples, the initial intervals were given as $[-10^8,10^8]$.

\begin{table}
\small
\centering
\begin{tabular}{||r|r|r|r||}
\hline
\hline
- 0.249150680 & + 0.125016350 & - 0.635550070  & + 1.48947730  \\
+ 1.609135400 & - 0.686607360 & - 0.115719920  & + 0.23062341  \\
+ 0.279423430 & - 0.119228120 & - 0.666404480  & + 1.32810730  \\
+ 1.434801600 & - 0.719940470 & + 0.110362110  & - 0.25864503  \\
+ 0.000000000 & - 0.432419270 & + 0.290702030  & + 1.16517200  \\
+ 0.400263840 & + 0.000000000 & + 1.258776700  & - 0.26908494  \\
- 0.800527680 & + 0.000000000 & - 0.629388360  & + 0.53816987  \\
+ 0.000000000 & - 0.864838550 & + 0.581404060  & + 0.58258598  \\
+ 0.074052388 & - 0.037157270 & + 0.195946620  & - 0.20816985  \\
- 0.083050031 & + 0.035436896 & - 1.228034200  & + 2.68683200  \\
- 0.386159610 & + 0.085383482 & + 0.000000000  & - 0.69910317  \\
- 0.755266030 & + 0.000000000 & - 0.079034221  & + 0.35744413  \\
+ 0.504201680 & - 0.039251967 & + 0.026387877  & + 1.24991170  \\
- 1.091628700 & + 0.000000000 & - 0.057131430  & + 1.46773600  \\
+ 0.000000000 & - 0.432419270 & - 1.162808100  & + 1.16517200  \\
+ 0.049207290 & + 0.000000000 & + 1.258776700  & + 1.07633970  \\
+ 0.049207290 & + 0.013873010 & + 2.162575000  & - 0.69686809  \\
\hline
\hline
\end{tabular}
\caption{Coefficients for the Inverse Kinematics Example.}
\label{kin-coef}
\end{table}

The results of {\tt Newton} on these two benchmarks are given in Table
\ref{kin-newton}. {\tt Newton} is fast on the first benchmark and
does not branch  much to obtain all solutions. The algorithm in
\cite{Hong95} branches more (the reported figure is 257 branches but
it is not really comparable due to the nature of the algorithm) and is
about 16 times slower on comparable machines. We are not aware of the
results of continuation methods on this problem. {\tt Newton} is
slower on the second application and takes about 6 minutes.  The
continuation method described in \cite{Cools94} requires about 30
seconds on a {\tt DEC 5000/200}. This method exploits the fact that
the Newton polytopes for the last 4 equations are the same. Note that
{\tt HRB} requires about 1630 and 4730 seconds on these examples. Note
also that {\tt Newton} can prove the existence of solutions in the
final intervals for these problems.

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r||}
\hline
\hline
               &{\tt kin1} &{\tt kin2} \\
\hline
{\it CPU time} &  14.240   &  353.06 \\
{\it branching}   &  89       &  5693  \\
\hline 
{\it na-ee}    &  17090    & 784687 \\
{\it na-te}    &  10176    & 123032 \\
{\it na-tot}   &  27266    & 907719 \\
\hline   
{\it fe-ee}    &  45656    & 1714779 \\
{\it fe-te}    &  62080    & 854384 \\
{\it fe-tot}   &  107736   & 2569163 \\
\hline 
{\it pr-con}   &  163      & 9505 \\
\hline 
{\it proof}    &  yes     & yes    \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Kinematics Benchmarks}
\label{kin-newton}
\end{table}



\subsection{An Economics Modelling Application}

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r|r||}
\hline
\hline
               &  4   &  5    &  6     &   7   &   8     &  9   \\
\hline
{\it CPU time} & 0.21 & 1.22  & 8.20   & 46.59 & 352.80  & 3311.42 \\
{\it growth}   &      & 5.80  & 6.72   & 5.68  & 7.57    & 9.38  \\
{\it branching}   & 24   & 119   & 517    & 2231  & 12248   & 82579 \\
\hline 
{\it na-ee}    & 834  & 4701  & 42481  & 214430& 1622417 & 14031838 \\
{\it na-te}    & 480  & 1976  & 6325   & 29238 & 157402  & 1219960 \\
{\it na-tot}   & 1314 & 6677  & 48806  & 243668& 1779819 & 15251798 \\
\hline   
{\it fe-ee}    & 2220 & 116628& 99994  & 489894& 3626847 & 30782836 \\
{\it fe-te}    & 1293 & 6304  & 29825  & 160284& 1062789 & 9118448 \\
{\it fe-tot}   & 3513 & 17932 & 129819 & 650178& 4689636 & 39901284 \\
\hline 
{\it pr-con}   & 37   & 147   & 687    & 2828  & 15265   & 104352 \\
\hline 
{\it proof}    & yes  & yes   & yes    & yes   & yes     & yes \\
\hline
\end{tabular}
\caption{{\tt Newton} on the economics modelling problem with initial intervals in $[-100,100]$}
\label{eco-newton}
\end{table}

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r|r|r||}
\hline
\hline
               &  4    &  5    &  6     &   7   &   8     &  9   \\
\hline
{\it CPU time} & 0.60  & 3.35  & 22.53  & 127.65 &  915.24 & 8600.28 \\
{\it growth}   &       & 5.58  & 6.72   & 5.66   &  7.16   & 9.39   \\    
{\it branching}   & 102   & 500   & 1778   & 7527   &  38638  & 244263 \\
\hline  
{\it na-ee}    & 2689  & 15227 & 122662 &  606805& 4413150 & 36325819 \\
{\it na-te}    & 978   & 3296  & 14055  & 64632  & 366436  & 2647016 \\
{\it na-tot}   & 3667  & 18523 & 136717 & 671437 & 4779586 & 38972835 \\
\hline    
{\it fe-ee}    & 7140  & 38160 & 291206 & 1400216& 9913850 & 80264897 \\
{\it fe-te}    & 3288  & 15024 & 80110  & 429396 & 2825368 & 22672208 \\
{\it fe-tot}   & 10428 & 53184 & 371716 &1829612 & 12739218& 102937105 \\
\hline 
{\it pr-con}   & 148   & 527   & 2080   & 8337   & 42704   & 271534  \\
\hline 
{\it proof}    & yes   & yes   & yes    & yes    & yes     & yes \\
\hline
\end{tabular}
\caption{{\tt Newton} on the economics modelling problem with initial intervals in $[-10^8,10^8]$}
\label{eco-newton-W}
\end{table}

\noindent The following example is taken from \cite{Morgan87}. It is a
difficult economic modelling problem that can be scaled up to
arbitrary dimensions. For a given dimension $n$, the problem can be
stated as the system
\[
\left\{ \begin{array}{l}
(x_k + \sum_{i = 1}^{n-k-1} x_i x_{i+k}) x_n - c_k = 0 \;\;\; (1
\leq k \leq n-1) \\
\sum_{l=1}^{n-1} x_l + 1 = 0
\end{array} \right.
\]
and the constants can be chosen at random. 

Table \ref{eco-newton} reports the results for various values of $n$
with an initial interval of $[-100,100]$. It is interesting to compare
those results with the continuation methods presented in
\cite{Cools94}. \cite{Cools94} reports times (on a DEC-5000/200) of
about 1 second for $n = 4$, 6 seconds for $n=5$, 50 seconds for $n=6$
and 990 seconds for $n=7$. {\tt Newton} is substantially faster on
this problem than this continuation method, since it takes about 47
seconds for $n =7$. More importantly, the growth factor seems much
lower in {\tt Newton}. The continuation method has growths of about 8
and 20 when going from 5 to 6 and 6 to 7, while {\tt Newton} has
growths of about 6.72 and 5.68.  Table \ref{eco-newton-W} gives the
same results for initial intervals in $[-10^8,10^8]$. It is
interesting to note that the computation times increase by less than a
factor 3 and that the growth is consistent with the results on the
smaller intervals. Note also the existence proof of solutions for
these problems. 

\subsection{Combustion Application}

\noindent This problem is also from Morgan's book \cite{Morgan87} and
represents a combustion problem for a temparature of $3000^{\circ}$.
The problem is described by the following sparse systems of equations
\[        
\left\{ \begin{array}{l}
x_2 + 2 \ x_6 + x_9 + 2\ x_{10} = 10^{-5}  \\
x_3 + x_8 = 3 \ 10^{-5}  \\
x_1 + x_3 + 2 \ x_5 + 2 \ x_8 + x_9 + x_{10} = 5 \ 10^{-5}  \\
x_4 + 2 \ x_7 = 10^{-5}  \\
0.5140437 \ 10^{-7} \ x_5 = x_1^2  \\
0.1006932 \ 10^{-6} \ x_6  = 2 \ x_2^2  \\
0.7816278 \ 10^{-15} \ x_7 = x_4^2  \\
0.1496236 \ 10^{-6} \ x_8  = x_1\ x_3  \\
0.6194411 \ 10^{-7} \ x_9 = x_1\ x_2  \\ 
0.2089296 \ 10^{-14} \ x_{10} = x_1\ x_2^2 \\
\end{array} \right.
\]
which is typical of chemical equilibrium systems. Table
\ref{combustion-newton} describes the results of {\tt Newton} on for
the initial intervals $[-1,1]$ and $[-10^8,10^8]$. {\tt Newton}
behaves  well on this example, since the continuation method of
\cite{Cools94} takes about 57 seconds. Note once again that a
substantial increase in the size of the initial intervals only induces
a slowdown of about 2.5 for {\tt Newton}. Note also the proof of existence of
the solutions. 

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r||}
\hline
\hline
               &  [-1,1]  & [-$10^8$,$10^8$] \\
\hline
{\it CPU time} &  4.06    & 9.94  \\
{\it branching}   &  183     & 523   \\
\hline  
{\it na-ee}    & 10676    & 28473 \\
{\it na-te}    & 5248     & 7952  \\
{\it na-tot}   &15924     & 36425 \\
\hline    
{\it fe-ee}    & 28928    & 77508 \\
{\it fe-te}    & 22464    & 49632 \\
{\it fe-tot}   & 51392    & 127140 \\
\hline 
{\it pr-con}   & 187      & 527 \\
\hline 
{\it proof}    & yes      & yes   \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Combustion Problem.}
\label{combustion-newton}
\end{table}

\subsection{Chemical Equilibrium Application}

\begin{table}
\small
\centering
\begin{tabular}{||r|r||r|r|r||r|r|r||r|r||}
\hline
\hline
{\it CPU time} & {\it branching} & {\it na-ee} & {\it na-te} & {\it na-tot} &
{\it fe-ee}    & {\it fe-te} & {\it fe-tot} & {\it pr-con}  & {\it
  proof} \\
\hline
6.32 & 256 & 13725 & 3400 & 17125 & 34811 & 17425 & 52236 & 425 & yes \\
\hline
\hline
\end{tabular}
\caption{{\tt Newton} on the Chemestry Problem with Initial Intervals
  $[0,10^8]$}
\label{chemestry-newton}
\end{table}

\noindent This problem originates from \cite{Morgan90} and describes a
chemical equilibrium system. The set of equations is as follows:
\[
\left\{ \begin{array}{l}
R = 10  \\
R_5 = 0.193  \\
R_6 = 0.002597 / \sqrt{40}  \\
R_7 = 0.003448 / \sqrt{40} \\
R_8 = 0.00001799 / 40 \\
R_9 = 0.0002155 / \sqrt{40} \\
R_{10} = 0.00003846 / 40  \\
x_1 \ x_2 + x_1 - 3 \ x_5 = 0  \\
2 \ x_1 \ x_2 + x_1 + x_2 \ x_3^2 + R_8 \ x_2 - R \ x_5 + 2 \ R_{10} \ x_2^2 +  R_7 \ x_2 \ x_3 + R_9 \ x_2 \ x_4  = 0  \\
2 \ x_2 \ x_3^2 + 2 \ R_5 \ x_3^2 - 8 \ x_5 + R_6 \ x_3 + R_7 \ x_2 \ x_3 = 0  \\
R_9 \ x_2 \ x_4 + 2 \ x_4^2 - 4 \ R \ x_5 = 0  \\  
x_1 \ x_2 + x_1 + R_{10} \ x_2^2 + x_2 \ x_3^2 + R_8 \ x_2 + R_5 \ x_3^2 + 
x_4^2 - 1 + R_6 \ x_3 + R_7 \ x_2 \ x_3 + R_9 \ x_2 \ x_4 = 0
\end{array}
\right.
\]
and all $x_i$ must be positive. The results are depicted in Table
\ref{chemestry-newton} for an initial interval $[0,10^8]$.
They indicate that {\tt Newton} is particularly effective on this
problem, since it takes about 6 seconds and proves the existence of a
solution in the final intervals.  Note that the continuation method of
\cite{Cools94} takes about 56 seconds on this problem.


\subsection{A Neurophysiology Application}

We conclude the experimental section by showing an example
illustrating the limitations of {\tt Newton}. The application is from
neurophysiology \cite{Cools94} and consists of the following system of
equations:
\[
\left\{ \begin{array}{l}
x_1^2 + x_3^2 = 1  \\
x_2^2 + x_4^2 = 1  \\        
x_5  \ x_3^3 + x_6 \ x_4^3 =  c_1  \\
x_5 \ x_1^3 + x_6 \ x_2^3 = c_2  \\
x_5 \ x_1 \ x_3^2 + x_6 \ x_4^2 \ x_2 = c_3  \\ 
x_5 \ x_1^2 \ x_3 + x_6 \ x_2^2 \ x_4 = c_4
\end{array} \right.
\]
No initial intervals for the variables were given and the constants
$c_i$ can be chosen at random. The continuation method of \cite{Cools94}
solves this problem in about 6 seconds. The results of {\tt Newton}
are depicted in Table \ref{neuro-newton} for various initial
intervals. {\tt Newton} is fast when the initial intervals are small
(i.e., $[-10,10]$). Unfortunately, the running time of the algorithm
increases linearly with the size of the initial intervals, showing a
limitation of the method on this example.

\begin{table}
\small
\centering
\begin{tabular}{||l||r|r|r|r||}
\hline
\hline
               &  $[-10,10]$  & $[-10^2,10^2]$ & $[-10^3,10^3]$ &
               $[-10^4,10^4]$ \\
\hline
{\it CPU time} &  0.91  & 11.69  & 172.71  & 2007.51  \\
{\it growth}   &        & 12.84  & 14.77   & 11.62    \\
{\it branching}   &  52    & 663    & 9632    & 115377   \\
\hline  
{\it na-ee}    & 4290   & 57224  & 645951  & 6541038  \\
{\it na-te}    & 810    & 6708   & 96888   & 1173456  \\
{\it na-tot}   & 5100   & 63932  & 742839  & 7714494  \\
\hline    
{\it fe-ee}    & 11012  & 144843 & 1620538 & 16647442 \\
{\it fe-te}    & 4104   & 48804  & 769056  & 9376632  \\
{\it fe-tot}   & 15116  & 193647 & 2389594 & 26024074 \\
\hline 
{\it pr-con}   & 69     & 983    & 15980   & 195270   \\
\hline 
{\it proof}    & yes    & yes    & yes     & yes      \\
\hline
\end{tabular}
\caption{{\tt Newton} on the Neurophysiology Problem.}
\label{neuro-newton}
\end{table}

\section{Related Work and Discussion}
\label{related}

The research described in this paper originated in an attempt to
improve the efficiency of constraint logic programming languages based
on intervals such as {\tt BNR-Prolog} \cite{BNR} and {\tt
  CLP(BNR)} \cite{CLP-BNR}. These Constraint Logic Programming
languages use constraint solving as basic operation and they were
based on the simple generalization of arc-consistency described
previously, i.e.
\[
I_i = \bo{I_i \cap \{ \ r_i \ \mid \ \exists r_1 \in I_1, \ldots, \exists
r_{i-1} \in I_{i-1}, \ldots, \exists r_{i+1} \in I_{i+1}, \ldots,
\exists r_{n} \in I_n: c(r_1,\ldots,r_n) \ \}}.
\]
This approximation was enforced on simple constraints such as 
\[
x_1 = x_2 + x_3, \;\; x_1 = x_2 - x_3, \;\; x_1 = x_2 \times x_3
\]
and complex constraints were decomposed in terms of these simple
constraints. 

As mentioned previously, this approach is not very effective and our
main goal was to design new approximations of arc-consistency that
could make use of existing interval Newton methods. The main problem
was the difficulty in characterizing the pruning of the Newton
operator $N^*$ in a declarative way (in order to introduce it nicely
in the above programming languages) and box-consistency emerged as an
attempt to generalize the operator to make sure that the bounds of the
interval were locally consistent. Subsequent research made us realize
that box-consistency is independent of the Newton operator and can be
enforced even if the functions are not continuous or differentiable.
In addition, the value of applying box-consistency on several
extensions became clear. On the one hand, box-consistency on the
Taylor extension generalizes interval methods based on Gauss-Seidel
iterations and enables us to capture nicely Hansen-Segupta's operator.
On the other hand, box-consistency on the natural and expanded
extensions is really orthogonal to the pruning obtained from the
Taylor expansion, producing a particularly effective algorithm. It is
also worth pointing out that {\tt Newton} spends most time in the
natural and expanded extensions. However, for many applications, the
use of the Taylor interval extension is critical to terminate the
search quickly and to avoid generating many small intervals around the
solutions. As a general observation, box-consistency on the natural
and expanded extensions seem effective when far from a solution while
box-consistency on the Taylor expansion seems effective when near a
solution. It is worth mentioning that that the interval community
has spent much effort to design additional techniques to speed up
further the computation when near a solution but have not considered
techniques to improve pruning when far from a solution.

It is interesting to note that the idea of using approximations of
arc-consistency was also used independently by Hong and Stahl
\cite{Hong95}, who were also exposed to research on Constraint Logic
Programming. Their use of projections is however quite different from
ours. The key idea is to work with a set of boxes and to use
projections to split a box into several subboxes by isolating all
zeros of a projection. This gives an algorithm of a very different
nature which cannot easily be characterized as a branch \& prune
algorithm since constraints are used to branch. Our approach seems to
be more effective in practice, since their use of projections may
generate many subboxes that may all need to be pruned away later on,
implying much redundant work. Our approach postpones the branching
until no pruning takes place and generates only subboxes when they are
strictly necessary to progress. It is also very interesting to report
that on all benchmarks that we tested the projection never isolates
more than two zeros, showing that box-consistency may be the right
trade-off here.  Finally, note that their approach seems to use
implicitly an expanded extension\footnote{The idea of sandwitching the
  interval function in between two real functions is described there.}
but they do not make use of the natural extension which is very
important for some applications.

The research described here also provides a uniform framework to
integrate these techniques in Constraint Logic Programming, to
understand the importance of the various pruning operators and their
relationships and to suggest further research directions. For
instance, higher notions of consistency such as path-consistency
\cite{Montanari74} may be worth investigating for some applications.

\section{Conclusion}
\label{conclusion}

In this paper, we presented a branch \& prune algorithm to find all
isolated solutions to a system of polynomial constraints over the
reals. The algorithm is based on a single concept, box-consistency,
which is an approximation of arc-consistency, a notion well-known in
artificial intelligence. Box-consistency can be instantiated by to
produce Hansen-Segupta operator as well as other narrowing operators
which are more effective when the computation is far from a solution.
The algorithm and its mathematical foundations are simple. Moreover,
the algorithm is shown to behave well on a variety of benchmarks from
kinematics, mechanics, chemistry, combustion, and economics. It
outperforms the interval Newton methods we know of and compares well
with continuation methods on their benchmarks. In addition, problems
such as the Broyden banded function and the Mor\'{e}-Cosnard
discretization of a nonlinear integral equation can be solved for
several hundred variables.  Limitations of the method (e.g. a
sensitivity to the size of the initial intervals on some problems)
have also been identified.

\section*{Acknowledgments}

We would like to thank F. Benhamou, A. Colmerauer, and B. Le Charlier
for many interesting discussions on this topic.  This research was
partly supported by the Office of Naval Research under grant
N00014-91-J-4052 ARPA order 8225, the National Science Foundation
under grant numbers CCR-9357704, a NSF National Young Investigator
Award.

\bibliographystyle{plain}
\bibliography{pascal}

\newpage

\appendix


\end{document}

Note that working with unions of intervals is called extended interval
arithmetic in the interval arithmetic literature.






\begin{figure}
\begin{tt}
\begin{tabbing}
123\=123\=123\=123\=123\=123\=\kill
{\bf function} {\it ShrinkL}( $F$: $\I \rightarrow \I$; $F'$: $\I
\rightarrow \I$; $I: \I$)$\fto \I$; \\
{\bf begin} \\
\>      success := false; \\
\>      rightI := {\it right(I)}; \\
\>      {\it {\tt PUSH}(I)}; \\
\>      {\bf while}  $\neg$ EMPTY\_STACK $\wedge$ $\neg$ success {\bf do} \\
\>\>            {\it I} := POP\_STACK; \\
\>\>            $I$ := $N^*(F,F',I)$; \\
\>\>            {\bf if} $0 \in F(I)$ {\bf then} \\
\>\>\>              {\bf if} {\it LeftBounded(F,I)} {\bf then} \\
\>\>\>\>                   success := true \\
\>\>\>              {\bf else} \\
\>\>\>\>                   PUSH([{\it center(I)},{\it right(I)}]); \\
\>\>\>\>                   PUSH([{\it left(I)},{\it center(I)}]);\\
\>\>\>              {\bf endif} \\
\>\>            {\bf endif} \\
\>      {\bf endwhile} \\
\>      {\bf if} $\neg$ success {\bf then} {\bf return} $\emptyset$ {\bf else return} [{\it left(I)},rightI] \\
{\bf end}
\end{tabbing}
\end{tt}
\caption{Implementation of Function {\it ShrinkL} using Newton Iteration}
\label{shr}
\end{figure}



In the following, we use an operation $\UDIV(I_c,I_a,I_d)$ which
computes $I_c - I_a/I_d$ using unions of intervals. In addition, we
take the convention that $\UDIV(I_c,I_a,I_d)$  always returns
a union $I_n \cup I_p$, where $I_n$ is empty or contains only
nonpositive numbers and $I_p$ is empty or contains only nonnegative
numbers. 



\begin{definition}[Newton Reduction] The Newton
  reduction operator $N: ((\I \fto \I) \times (\I \fto \I) \times \I) \fto
  (\I \times \I)$ 
\[
N(F,F',I) = \UDIV(\ap{c},F(\ap{c}),F'(I)).
\]
where $c = \CENTER(I)$. The iterated Newton reduction operator 
$N^*: (\I \fto \I \times \I \fto \I \times \I) \fto \I$  is defined as
\begin{tt}
\begin{tabbing}
12345\= \kill
$N^*(F,F',I) = I_n$  $(n \geq 1)$ where \\
\>  $I_0 = I$ \\
\>  $I_n \cup I_p = \UDIV(F,F',I_i)$ \\
\>  $I_{i+1} = (I_i \cap I_n) \cup (I_i \cap I_p)$ \\
\>  $I_n = I_{n-1}$.
\end{tabbing}
\end{tt}
\end{definition}

The following property is a condition to guarantee the correctness of
the Newton reduction operators. 
\begin{definition}[Newton Compatibility] $F': \I \fto \I$ is
  are Newton-compatible with $F: \I \fto \I$ and $I$ iff
\[
\forall I' \subseteq I \ \forall r \in I': \ 0 \in F(\ap{r})
\Rightarrow r \in N(F,F',I').
\]
\end{definition}

\noindent Assuming a function $F'$ that is Newton-compatible with $F$
and $I$, the implementation of $ShrinkL$ can be refined as shown in
Figure \ref{shr}. The only change is the addition of the iterated
Newton reduction operator to reduce the interval which is removed from
the stack.


\begin{quote}
$[-0.3894150838568418705,-0.3894150838568416483]$ \\
$[-0.9210623716474195533,-0.9210623716474192201]$ \\
$[+0.5646048742322073987,+0.5646048742322289372]$ \\
$[-0.8253613366236840499,-0.8253613366236680626]$ \\
$[+0.7174034286244790337,+0.7174034286245180026]$ \\
$[-0.6966579652870119599,-0.6966579652869752115]$ \\
$[+0.8414785200543410947,+0.8414785200543567490]$ \\
$[-0.5402905702371235686,-0.5402905702370986995]$ \\
$[+0.9320048415873238845,+0.9320048415873262160]$ \\
$[-0.3624458238934301080,-0.3624458238934238907]$ \\
$[-0.8559151179237649565,-0.8559151123762789303]$ \\
$[-0.5171163509365193090,-0.5171163417544921703]$
\end{quote}

When $n=5$, {\tt Newton} returns the intervals
\begin{quote}
$[-0.4283028646427008979,-0.4283028646427007313]$ \\
$[-0.4765965315010954329,-0.4765965315010952663]$ \\
$[-0.5196377221007547620,-0.5196377221007544289]$ \\
$[-0.5588619565270254164,-0.5588619565270250832]$ \\
$[-0.5588619565270254164,-0.5588619565270250832]$
\end{quote}

