\documentstyle[12pt,amssymb]{amsart}

% \headheight=8pt     \topmargin=0pt
% \textheight=624pt   \textwidth=432pt
% \oddsidemargin=18pt \evensidemargin=18pt

% other margins :

 \headheight=8pt     \topmargin=-15pt
 \textheight=634pt   \textwidth=470pt
 \oddsidemargin=2pt  \evensidemargin=2pt   

% using PicTeX :

\input{prepictex}
\input{pictex}
\input{postpictex}

% some environments :

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{proposition}{Proposition}[section] 
\newtheorem{algorithm}{Algorithm}[section] 

% using AMS mathematical symbols

\newcommand{\nn}{{\Bbb N}}
\newcommand{\ee}{{\Bbb E}}
\newcommand{\zz}{{\Bbb Z}}
\newcommand{\rr}{{\Bbb R}}
\newcommand{\cc}{{\Bbb C}}
\newcommand{\pp}{{\Bbb P}}                 

\def\qed{\hfill $\Box$}  % end of proof  
\def\qex{\hfill \quad\vrule height1.2ex width0.5em depth 0pt} % end of example 

\begin{document}

\title{Toric Newton Method for Polynomial Homotopies}

\author{Jan Verschelde}
\address{\hskip-\parindent
		Mathematical Sciences Research Institute \\
        1000 Centennial Drive \\
		Berkeley, CA 94720-5070, USA.}
\address{\hskip-\parindent
        Department of Mathematics \\
        Michigan State University \\
        East Lansing, MI 48824-1027, USA.}

\email{jan@msri.org, jan@math.msu.edu or jan.verschelde@na-net.ornl.gov,
and urls: http://www.msri.org/people/members/jan 
http://www.mth.msu.edu/{\~{}}jan}

\thanks{Research at MSRI is supported by a post-doctoral fellowship.
Supported in part by NSF under Grant DMS - 9804846 at Michigan State 
University}

\bibliographystyle{plain}

\begin{abstract}
\noindent This paper defines a generalization of Newton's method to deal
with solution paths defined by polynomial homotopies that lead to extremal
values.  Embedding the solutions in a toric variety leads to explicit scaling
relations between coefficients and solutions.  Toric Newton is a 
symbolic-numeric algorithm where the symbolic pre-processing exploits the
polyhedral structures.  The numerical stage uses the additional variables
introduced by the homogenization to scale the components of the solution
vectors to the complex unit circle.  Toric Newton generates appropriate
affine charts and enables to approximate the magnitude of large solutions
of polynomial systems.

\medskip

\noindent {\bf AMS Subject Classification.} 14M25, 52B55, 65H10, 68Q40.
% 14M25 Toric varieties
% 52B55 Computational aspects related to convexity
% 65H10 Systems of nonlinear equations
% 68Q40 Symbolic Computation, algebraic computation

\medskip

\noindent {\bf Keywords}.  polynomial system, toric varieties,
Newton's method, homotopy continuation.

\end{abstract}

\maketitle

\section{Introduction}

If a polynomial system has a large root, we first want to estimate
its magnitude before the approximation of its actual value.
In the context of polyhedral homotopies we can certify~\cite{hv98} diverging
solution paths and separate those from paths leading to extremal values.
The purpose of this paper is to show how polyhedral methods can approximate
the magnitude of large roots.
In~\cite{hs95} and~\cite{vvc94} polyhedral homotopies were introduced
to exploit sparsity.  The number of solution paths equals the bound of
Bernshte\v{\i}n's theorem.
First we give some notations to state Bernshte\v{\i}n's theorems.

\smallskip

Denote ${\bf x}^{\bf a} = x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n}$.
The {\em support} of a polynomial 
$f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}$ is the
set $A = \{ \ {\bf a} \in \zz^n \ | \ c_{\bf a} \not= 0 \ \}$.
The {\em Newton polytope} $P$ of $f$ is the convex hull of its support~$A$.
The structure of a polynomial system $F({\bf x}) = {\bf 0}$ is represented
by the Newton polytopes ${\cal P} = (P_1,P_2,\ldots,P_n)$ spanned by
the supports ${\cal A} = (A_1,A_2,\ldots,A_n)$.
Denote $\cc^* = \cc \setminus \{ 0 \}$.
In~\cite{ber75} Bernshte\v{\i}n proved that the number of isolated solutions
in~$(\cc^*)^n$ of a polynomial system with generic coefficients 
equals the mixed volume $V_n({\cal P})$ of its Newton polytopes~$\cal P$.

\smallskip

In the nongeneric case we examine faces of systems.
A {\em face}~$\partial_{\bf v} P$ of a polytope~$P$ is characterized by a 
direction~${\bf v} \in \rr^n \setminus \{ {\bf 0} \}$ and spanned by
$\partial_{\bf v} A = \{ \ {\bf a} \in A \ | \ \forall {\bf b} \in A:
\langle {\bf a} , {\bf v} \rangle \leq \langle {\bf b} , {\bf v} \rangle \ \}$.
A {\em face polynomial of~$f$} is~$\partial_{\bf v} f({\bf x}) =
 \sum_{{\bf a} \in \partial_{\bf v} A} c_{\bf a} {\bf x}^{\bf a}$.
A {\em face system of~$F$} is $\partial_{\bf v} F({\bf x}) = {\bf 0}$ with
equations $\partial_{\bf v} f_i({\bf x}) = 0$, $i=1,2,\ldots,n$.
For~$V_n({\cal P}) > 0$,
Bernshte\v{\i}n's second theorem~\cite{ber75} states that
$F({\bf x}) = {\bf 0}$ has fewer than~$V_n({\cal P})$ isolated solutions
in~$(\cc^*)^n$ if and only if there is some
$\partial_{\bf v} F({\bf x}) = {\bf 0}$ that has a solution in~$(\cc^*)^n$.

\smallskip

A toric variety $X$ is a solution set of a polynomial system that contains
the torus $(\cc^*)^n$ as a dense open subset with $(\cc^*)^N$ as a group
acting on~$X$.
In section~2 we make this definition from~\cite{ful93} concrete in
applying the homogenization procedure described in~\cite{cls98}.
The group action is used in the third section to scale the solution
components to the complex unit circle, narrowing the torus $(\cc^*)^n$
to a numerically favorable region.
See~\cite{cox95}, \cite{dan78}, \cite{ful93}, \cite{gkz94}, and~\cite{oda88}
for more on toric varieties.

\smallskip

Representing diverging and badly-scaled solutions is a common struggle
in the daily life of a numerical analysist.
In~\cite[Chapter 5]{mor87}, \cite{mm87}, \cite{wbm87,wsmmw97}
we find coefficient and equation scaling,
combined with centering, an effective and useful tool.
The classical homogenization~\cite{bcss97}, \cite{li97}, \cite{mor87}
in projective space easily leads to clustering and fails to recognize
the deficiency behavior when applied to sparse polynomial systems.
In the third section of this paper we exemplify this situation.

\smallskip

In the fourth section we define Toric Newton method and analyze the cost
of its execution.  This algorithm is a symbolic-numeric approach.
The Toric Newton method consists of a symbolic pre-processing stage,
where the homogeneous forms are determined and the normal fans are
refined to be simplicial.  This stage is performed only once.
The numerical stage is executed iteratively while following the solution
paths and takes place in the original space of dimension~$n$.

\smallskip

Toric Newton augments the robustness of polyhedral end games~\cite{hv98}
as used in numerical homotopies~\cite{li97}.
Convergence properties of Toric Newton are addressed in the fifth section.
The complexity of Newton in the dense case has been analyzed in~\cite{bcss97}
and has recently been extended to sparse complexity models~\cite{ded97},
multi-homogeneous~\cite{ds97} and over-determined systems~\cite{ds98}.
The Toric Newton method parallels similar developments in providing
algorithms using torus embeddings to compute sparse resultants~\cite{roj98},
residues~\cite{cd97},~\cite{cds96}, and Bezoutians~\cite{rs98}.
All these methods rely on polytopes~\cite{zie95}.
Here we give practical evidence for the performance of the polyhedral
constructions.

\section{Newton Polytopes and Homogeneous Coordinates}

If the Newton polytope lies in a hyperplane, then the coordinates of
the polynomial are homogeneous.  The exponents of the scaler $\lambda$ 
are given by the components of the direction that is perpendicular to
the hyperplane, as expressed formally in the following proposition.

\begin{proposition} {\rm (\cite[page~193]{gkz94})} Consider
${\displaystyle f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}}$.
\newline If $\exists {\bf v} \in \zz^n \setminus \{ {\bf 0} \}$,
$\exists d \in \zz$:
$\forall {\bf a} \in A, \langle {\bf a} , {\bf v} \rangle = d$,
then $f( \lambda^{v_1} x_1 , \lambda^{v_2} x_2 , \ldots , \lambda^{v_n} x_n )
= \lambda^d f({\bf x})$.
\end{proposition}
{\em Proof.}  Denote $\lambda^{\bf v} {\bf x} = \lambda^{v_1} x_1 
\lambda^{v_2} x_2 \cdots \lambda^{v_n} x_n$,
then $(\lambda^{\bf v} {\bf x})^{\bf a} = (\lambda^{v_1} x_1)^{a_1}
 (\lambda^{v_2} x_2)^{a_2} \cdots ( \lambda^{v_n} x_n)^{a_n}
  = \lambda^{\langle {\bf a} , {\bf v} \rangle}$. 
So, ${\displaystyle f(\lambda^{\bf v} {\bf x})
 = \sum_{{\bf a} \in A} c_{\bf a} (\lambda^{\bf v} {\bf x})^{\bf a}
 = \sum_{{\bf a} \in A} c_{\bf a} \lambda^{\langle {\bf a}, {\bf v} \rangle}
                                                      {\bf x}^{\bf a}
 = \lambda^d f({\bf x})}$.~\qed

\medskip

The multilinearity of the inner product provides us with a whole range
of homogeneous coordinates when the polytope lies in the
intersection of many hyperplanes.

\begin{example} \label{exproj} (projective space)
{\rm Let $(x,y,z) \in \pp^2$ and
consider $f(x,y,z) = x+y+z$.  The support of $f$ is the set
$A = \{ (1,0,0) , (0,1,0) , (0,0,1) \}$ and lies in the hyperplane
$\langle {\bf a} , {\bf v} \rangle = d$ with ${\bf v} = (1,1,1)$, $d = 1$.
The homogeneous coordinates 
are $(\lambda^1 x , \lambda^1 y , \lambda^1 z )$.~\qex }
\end{example}

\begin{example} \label{exmulti} (multi-projective space) {\rm
Let $(x_0,x_1,y_0,y_1) \in \pp^1 \times \pp^1$ and consider
$f(x_0,x_1$, $y_0,y_1) = x_1 y_1 - x_0 y_0$.  The support of $f$ is the set
$A = \{ (0,1,0,1) , (1,0,1,0) \}$ and lies in the hyperplanes
$\langle {\bf a} , {\bf v}^i \rangle = d$ with ${\bf v}^1 = (1,1,0,0)$,
${\bf v}^2 = (0,0,1,1)$, ${\bf v}^3 = (1,0,0,1)$, and $d = 1$.
First we have the two usual bi-homogeneous coordinates
$(\lambda^1 x_0 , \lambda^1 x_1 , y_0,y_1)$
and $(x_0,x_1,\lambda^1 y_0 , \lambda^1 y_1 )$, which allow to scale
the $x$-variables independently from the $y$'s.
In addition, we have the third relation, because the Newton polytope is 
an edge and not of dimension two as one would expect.
The corresponding homogeneous coordinates are
$(\lambda^1 x_0, x_1 , y_0 , \lambda^1 y_1)$.~\qex }
\end{example}

Examples~\ref{exproj} and~\ref{exmulti} are particular cases of
homogeneous coordinates.  The general definition is given below.
In applying the homogenization method in~\cite{cls98} we copy
the original monomials~${\bf x}^{\bf a}$.
By this adaption we recognize the structure of the polytope and have
explicit scaling relations between the coefficients and solutions,
as needed in the toric version of Newton's method.

\begin{definition} \label{deftorhom} {\rm Consider
${\displaystyle f({\bf x}) = \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a}}$,
with $P = {\rm conv}(A)$.  Assume ${\rm vol}(P) > 0$.
Let $P$ have $N$ facets given by $N$ inward pointing normals
${\bf v}^i \in \zz^n \setminus \{ {\bf 0} \}$ with all components
relatively prime.  For every ${\bf v}^i$,
define $d_i = \min_{{\bf a} \in A}  \langle {\bf a} , {\bf v}^i \rangle$.
The facets are spanned by
$\partial_{{\bf v}^i} A = \{ {\bf a} \in A | 
  \langle {\bf a} , {\bf v}^i \rangle = d_i \}$.
The {\em homogenization map $\phi_A$} introduces for every ${\bf v}^i$ 
a new variable $z_i$, as follows:
\begin{equation} \label{eqdeftorhom}
  \phi_A ( {\bf x}^{\bf a} ) = {\bf x}^{\bf a}
    \prod_{i=1}^N z_i^{\langle {\bf a} , {\bf v}^i \rangle - d_i}
\quad {\rm and} \quad
  \phi_A \left( \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a} \right)
  = \sum_{{\bf a} \in A} c_{\bf a} \phi_A ( {\bf x}^{\bf a} ).
\end{equation}
If ${\rm vol}(P) = 0$, then $A$ lies in some hyperplanes
$\langle {\bf x} , {\bf v}^i \rangle = c_i$, and we add $\pm {\bf v}^i$
to the normals of the $k$-faces of $P$, where $k = {\rm dim}(A) - 1$.  }
\end{definition}
The result of the homogenization is the map
$\phi_A(f) : (\cc^*)^n \times \cc^N \rightarrow \cc
           : ({\bf x}, {\bf z}) \mapsto \phi_A(f)({\bf x}, {\bf z})$.
We may allow zero values for the $z$-variables, because of the
following property.

\begin{proposition} \label{prophomcoord1}
Let $\deg(f,x)$ denote the maximal exponent of the variable $x$
in the polynomial $f$.  Then,
\begin{equation}
\deg\left(\phi_A \left(
   \sum_{{\bf a} \in A} c_{\bf a} {\bf x}^{\bf a} \right) , z_i \right) \geq 0
\quad {\rm and} \quad
\deg(\phi_A ( {\bf x}^{\bf a} ) , z_i ) 
\left\{ \begin{array}{lll}
           = & 0 & {\rm if} \ {\bf a} \in \partial_{{\bf v}^i} A \\
           > & 0 & {\rm if} \ {\bf a} \notin \partial_{{\bf v}^i} A
        \end{array}
\right. .
\end{equation}
\end{proposition}
{\em Proof.}  The directions perpendicular to the facets are given
by inward pointing normals, i.e.: for a facet
normal ${\bf v}^i$: $d_i \leq \langle {\bf a} , {\bf v}^i \rangle$,
$\forall {\bf a} \in A$.  Therefore all exponents of $z_i$
of~$\phi_A$ are nonnegative.~\qed

\medskip

This homogenization corresponds geometrically to assigning as exponents
of $z_i$ the lattice distances relative to the $i$th facet of~$P$.
If we collect for every monomial the exponents of the new $z$-variables
and store these vectors rowwise in a matrix, then we obtain the
vertex-facet incident matrix~\cite{zie95} of the polytope~$P$.
The normal fan $\Sigma_P$ collects all cones that are normal to the faces
of~$P$.  The generators of those normal cones are the normals to the facets.

\begin{example} \label{runningex} (running example)
{\rm Consider $f(x_1,x_2) = c_{10} x_1 + c_{01} x_2
+ c_{20} x_1^2 + c_{12} x_1 x_2^2$.
The Newton polytope has four facets, with inward pointing normals
given by ${\bf v}^1 = (1,1)$, ${\bf v}^2 = (0,1)$, ${\bf v}^3 = (-2,-1)$, and  
${\bf v}^4 = (1,-1)$.  In Figure~\ref{figlatt}, the
homogenizing along~${\bf v}^4$ is illustrated.

\begin{figure}[hbt]
\begin{center}
\begin{picture}(360,120)(0,0)

\put(0,0){
\begin{picture}(120,120)(0,0)

\put(30,20){\vector(0,1){80}}
\put(20,30){\vector(1,0){120}}

\put(30,60){\circle*{3}}
\put(60,30){\circle*{3}}
\put(90,30){\circle*{3}}
\put(60,90){\circle*{3}}

\thicklines
\put(30,60){\line(1,-1){30}}  \put(32,77){$z_4$}
\put(30,60){\line(1,1){30}}   \put(37,37){$z_1$}
\put(60,30){\line(1,0){30}}   \put(68,20){$z_2$}
\put(60,90){\line(1,-2){30}}  \put(78,60){$z_3$}

\put(45,75){\vector(1,-1){15}}

\setplotarea x from 0 to 60, y from 0 to 60
\setdashes <4pt>
\setlinear
\plot -20 10   80 110 /
\plot  10 10  100 100 /
\plot  40 10  120  90 /
\plot  70 10  140  80 /

\put(-23,4){${}_0$}
\put(7,4){${}_1$}
\put(37,4){${}_2$}
\put(67,4){${}_3$}

\put(15,60){$P$}

\end{picture}
}

\put(200,0){
\begin{picture}(120,120)(0,0)

\put(60,60){\circle*{3}}
\put(60,60){\vector(1,1){40}}   \put(102,102){${\bf v}^1$}
\put(60,60){\vector(0,1){50}}   \put(62,112){${\bf v}^2$}
\put(60,60){\vector(-2,-1){50}} \put(0,35){${\bf v}^3$}
\put(60,60){\vector(1,-1){50}}  \put(115,8){${\bf v}^4$}

\setplotarea x from 0 to 60, y from 0 to 60
\setdashes <4pt>
\setlinear
\plot  60 68   65 65 /           % cone between v1 and v2
\plot  60 73   70 70 /
\plot  60 80   75 75 /
\plot  60 87   80 80 /
\plot  60 94   85 85 /
\plot  65 65   68 52 /           % cone between v1 and v4
\plot  70 70   76 44 /
\plot  75 75   84 36 /
\plot  80 80   92 28 /
\plot  85 85  100 20 /
\plot  52 56   60 68 /           % cone between v3 and v2
\plot  44 52   60 73 /
\plot  36 48   60 80 /
\plot  28 44   60 87 /
\plot  20 40   60 94 /
\plot  52 56   68 52 /           % cone between v3 and v4
\plot  44 52   76 44 /
\plot  36 48   84 36 /
\plot  28 44   92 28 /
\plot  20 40  100 20 /

\put(105,65){$\Sigma_P$}

\end{picture}
}

\end{picture}

\caption{On the left we see that the exponents of $z_4$ are lattice distances
relative to the fourth facet of~$P$ with normal~${\bf v}^4$.
On the right the normal fan $\Sigma_P$ of $P$ is displayed.}
\label{figlatt}
\end{center}
\end{figure}

The homogenization is
$\phi_A(f)({\bf x},{\bf z}) = c_{10} z_3^2 z_4^2 x_1 + c_{01} z_2 z_3^2 x_2
+ c_{20} z_1 z_4^3 x_1^2 + c_{12} z_1^2 z_2^2 x_1 x_2^2$.~\qex }
\end{example}

The homogeneous coordinates correspond to the inward pointing
normals to the facets.

\begin{proposition} \label{prophomcoord2}
Consider the facet normal ${\bf v}^i$ with corresponding $d_i$.
\newline Denote $\lambda^{{\bf v}^i} {\bf x} {\bf z}
 = (\lambda^{v_1^i} x_1 , \lambda^{v_2^i} x_2 , \ldots , \lambda^{v_n^i} x_n ,
	z_1,\ldots,z_{i-1}, \lambda^{-1} z_i, z_{i+1}, \ldots, z_n)$.
\newline Then
$\phi_A(f)(\lambda^{{\bf v}^i} {\bf x} {\bf z})
   = \lambda^{d_i} \phi_A(f)({\bf x},{\bf z})$.
\end{proposition}
{\em Proof.}  First consider one monomial~${\bf x}^{\bf a}$:
\begin{eqnarray}
  \phi_A({\bf x}^{\bf a}) \left.
     \right|_{({\bf x},{\bf z}) = \lambda^{{\bf v}^i} {\bf x} {\bf z}}
 & = & {\bf x}^{\bf a} 
     \prod_{j=1}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} \left.
     \right|_{({\bf x},{\bf z}) = \lambda^{{\bf v}^i} {\bf x} {\bf z}} \\
 & = & (\lambda^{v_1^i} x_1 \lambda^{v_2^i} x_2 \cdots 
	 \lambda^{v_n^i} x_n)^{\bf a} 
     \left(
       \prod_{j=1,\not=i}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} 
     \right)
     (\lambda^{-1} z_i)^{\langle {\bf a} , {\bf v}^i \rangle - d_i} \\
 & = & \lambda^{ \langle {\bf a} , {\bf v}^i \rangle } {\bf x}^{\bf a}
       \lambda^{ - \langle {\bf a} , {\bf v}^i \rangle + d_i }
       \prod_{j=1}^N z_j^{\langle {\bf a} , {\bf v}^j \rangle - d_j} \\
 & = & \lambda^{d_i} \phi_A({\bf x}^{\bf a})
\end{eqnarray}
The multilinearity of $\phi_A$ completes the proof.~\qed

\medskip

The transition from toric to affine space is immediate by taking
$\lambda = z_i$ in every $\lambda^{{\bf v}^i} {\bf x} {\bf z}$.

\begin{example} (Example~{\rm \ref{runningex}} continued.) {\rm
The homogeneous coordinates corresponding to the facets are given by
\begin{equation} \label{eqhomcoord}
  \begin{tabular}{ccccccccccccc}
    ( & $\lambda^1 x_1$    & , & $\lambda^1 x_2$ & , 
		   & $\lambda^{-1} z_1$ & , & $z_2$ & , & $z_3$ & , & $z_4$ & ) \\
    ( & $\lambda^0 x_1$    & , & $\lambda^1 x_2$ & , 
		   & $z_1$ & , & $\lambda^{-1} z_2$ & , & $z_3$ & , & $z_4$ & ) \\
    ( & $\lambda^{-2} x_1$ & , & $\lambda^{-1} x_2$ & , 
		   & $z_1$ & , & $z_2$ & , & $\lambda^{-1}z_3$ & , & $z_4$ & ) \\
    ( & $\lambda^1 x_1$    & , & $\lambda^{-1} x_2$ & ,
           & $z_1$ & , & $z_2$ & , & $z_3$ & , & $\lambda^{-1}z_4$ & ) \\
  \end{tabular}
\end{equation}
for the respective normals ${\bf v}^i$, $i=1,2,3,4$.
Observe that as $\lambda \rightarrow 0$, the $x$-values in affine
space go either to zero or infinity.
We have a multitude of homogeneous coordinates by taking combinations,
see $\Sigma_P$ in Figure~\ref{figlatt}.
Positive combinations of any two vectors normal to adjacent
facets correspond to constructing a vector inside the cone 
$\sigma \in \Sigma_P$ that has these vectors among its generators.
The relations~(\ref{eqhomcoord}) illustrate how the group~$(\cc^*)^N$ 
acts on the variety.~\qex }
\end{example}

For general polynomial systems not all supports are equal.
Suppose there are $r$ different supports which we group into
${\cal A} = (A_1,A_2,\ldots,A_r)$ spanning the
Newton polytopes ${\cal P} = (P_1,P_2,\ldots,P_r)$.
The idea is to embed the system using $r$ distinct sets of homogenization
variables, which provides a more cost-effective generalization than taking
the Minkowski sum of the polytopes in~$\cal P$.
This idea works because any normal to the Minkowski sum is a sum of normals
to the components in the sum, whence we can capture all faces in this way.
Definition~\ref{deftorhom} is extended to systems with different supports
as follows, taking $r=n$ for notational convenience. 

\begin{definition}  \label{defmixtorhom} {\rm 
Let $\cal A$ be the supports of $F({\bf x}) = {\bf 0}$.
Suppose the polytope spanned by $A_i$ has $N_i$ facet normals.
Abbreviate ${\bf z} = ({\bf z}^1, {\bf z}^2, \ldots, {\bf z}^n)$.
The {\em mixed homogenization map $\phi_{\cal A}$} is
\begin{equation}
   \phi_{\cal A}(F)
     : (\cc^*)^n \times \cc^{N_1} \times \cc^{N_2} \times
         \cdots \times \cc^{N_n} \rightarrow \cc^n
     : ({\bf x}, {\bf z}) \mapsto \phi_{\cal A}(F)({\bf x}, {\bf z}),
\end{equation}
with $\phi_{\cal A}(F)({\bf x}, {\bf z}) = (\phi_{A_1}(f_1)({\bf x},{\bf z}^1),
\phi_{A_2}(f_2)({\bf x},{\bf z}^2)$, $\ldots$, 
$\phi_{A_n}(f_n)({\bf x},{\bf z}^n))$
and $\phi_{A_i}(f_i)$ as in~{\rm (\ref{eqdeftorhom})}.
For a homotopy $H({\bf x},t) = {\bf 0}$, the continuation parameter $t$
is treated as a constant in setting up the homogenized
homotopy~$\phi_{\cal A}(H)({\bf x},{\bf z},t) = {\bf 0}$. }
\end{definition}

For several polytopes sharing the same facet normal, we can use the
same corresponding $z$-variable.  With this economical use of $z$-variables,
Definition~\ref{defmixtorhom} yields the same homogenization map $\phi_A$ as
in Definition~\ref{deftorhom} when applied to ${\cal A} = (A,A,\ldots,A)$.
Since the general homogenization is also determined by inward pointing
normals ${\bf v}^i$, Propositions~\ref{prophomcoord1} and~\ref{prophomcoord2}
have a direct generalization to the case of different supports and
are omitted.

\begin{proposition}  \label{propfacesols1}
Suppose $A$ is the support of $f({\bf x})$.
\newline Let $({\bf x}^*,{\bf z}^*) \in (\cc^*)^n \times \cc^N$ $:$
$\phi_A(f)({\bf x}^*,{\bf z}^*) = 0$.
\begin{enumerate}
  \item If $z^*_i = 0$ and $z^*_j \not= 0$, $\forall j \not= i$,
        then $\partial_{{\bf v}^i}f({\bf x}) = 0$
        has a solution in~$(\cc^*)^n$.
  \item Let $I = \{ \ i \ | \ z^*_i = 0,
            \# (\cap \partial_{{\bf v}^i} A) \geq 2 \ \}$,
        $z^*_j \not= 0$, $\forall j \not\in I$, and
		${\bf v} = {\displaystyle \sum_{i \in I} {\bf v}^i}$.
\newline Then $\partial_{\bf v}f({\bf x}) = 0$ has a solution in~$(\cc^*)^n$.
\end{enumerate}
The solutions ${\bf x}^*$ of~$\partial_{\bf v} f$
correspond to solutions of~$f$ with components equal to~$0$ and/or~$\infty$.
\end{proposition}
{\em Proof}.  Although the second statement generalizes the first,
we distinguish two cases:
\begin{enumerate}
  \item Suppose there is only one $z^*_i = 0$ with corresponding
        facet normal~${\bf v}^i$.
Then all monomials in $\phi_A(f)$ whose exponent
vector does not lie on the $i$th facet vanish when $z_i = z^*_i = 0$.
For $z^*_j \not= 0$, $\forall j \not= i$, the homogeneous coordinates
$\lambda^{{\bf v}^j} {\bf x}^*{\bf z}^*$ (as in Proposition~\ref{prophomcoord2})
are used successively, each time setting $\lambda = z^*_j$,
$\forall j \not= i$, to obtain a solution in~$(\cc^*)^n$ of
$\partial_{{\bf v}^i} f({\bf x}) = {\bf 0}$.

   \item When several $z^*_i$'s are zero, more monomials in the system vanish.
The condition on $I$ guarantees that there are still enough
monomials left in~$\partial_{\bf v} f$ to admit a solution in $(\cc^*)^n$.
Successive application of Proposition~\ref{prophomcoord2} setting
$\lambda = z^*_j$, $\forall j \notin I$, yields a solution
in $(\cc^*)^n$ of $\partial_{\bf v} f({\bf x}) = {\bf 0}$.

\end{enumerate}
The homogeneous coordinates contain as exponents
for $\lambda$ the components of the normals. 
Any normal $\bf v$ has at least one component different from zero, 
call it $v_k$.  Setting $\lambda = z_i \rightarrow 0$ to transform to
affine space will lead to $x_k \rightarrow \infty$, if $v_k < 0$,
or $x_k \rightarrow 0$, in case $v_k > 0$.~\qed

\medskip

In~\cite{cls98},
solutions that have no $z$-variables equal to zero, or only in cases
covered by Proposition~\ref{propfacesols1} are called~{\em nontrivial}.
Indeed, without the conditions on~$I$ in Proposition~\ref{propfacesols1}, all
monomials could vanish.
The generalization of Proposition~\ref{propfacesols1} to polynomial systems
is immediate.

\begin{proposition}  \label{propfacesols}
Let $F({\bf x}) = {\bf 0}$ have supports~$\cal A$.
The nontrivial solutions of $\phi_{\cal A}(F)({\bf x},{\bf z}) = {\bf 0}$
with some $z$-variables equal to zero correspond to solutions of face systems.
\end{proposition}    

We close this section with some remarks on systems with different supports.
In the homogenization we care for all coefficients in the system,
therefore we distinguish between supports, rather than between polytopes.
Secondly, for sparse polynomials it often happens that the supports lie
in lower-dimensional planes.  In that case we treat the normal vectors that
define the hyperplanes that contain the supports as facet normals in addition
to the normal to the faces that define the lower-dimensional polytopes.
Finally, aside from this complication, working with different supports does
not alter the homogenization method.

\section{Scaling Coefficients and Solutions}

The variables introduced by the homogenization allow to scale 
$\phi_A(f)({\bf x},{\bf z}) = {\bf 0}$, so that every component of~$\bf x$
lies on the unit circle in~$\cc^*$, while keeping the $z$-variables bounded.

\begin{proposition} \label{propscaling}
Suppose ${\bf x} \in (\cc^*)^n$
satisfies~$\phi_A(f)({\bf x},{\bf z}) = {\bf 0}$
for fixed $\bf z$ with $0 < z_j \leq 1$, $\forall j = 1,2,\ldots,N$.
Let ${\bf v} = (-\log|x_1|,-\log|x_2|,\ldots,-\log|x_n|)$, with
$\log$ of base 10.  Then
\begin{enumerate}
  \item $\exists {\bf c} \in \rr_+^n$,
        $\exists J = (j_1,j_2,\ldots,j_n)$, $j_i \in \{ 1,2,\ldots,N \}:$
        ${\displaystyle {\bf v} = \sum_{i=1}^n c_i {\bf v}^{j_i}}$.
  \item The homogeneous coordinates
\begin{equation}
   \left(
     \left. \lambda^{\left( - \sum_{i=1}^n c_i {\bf v}^{j_i} \right)_k} x_k
     \right|_{k=1}^n, \ldots,
     \left. \lambda^{c_i} z_{j_i} \right|_{i=1}^n , \ldots
   \right)
\end{equation}
with $\lambda = 0.1$ scale $({\bf x},{\bf z})$ to
$({\widetilde {\bf x}},{\widetilde {\bf z}})$ with 
$|{\widetilde x}_k| = 1$ and $0 < {\widetilde z}_j \leq 1$.
\end{enumerate}
\end{proposition}
{\em Proof}.  The above statements are true because:
\begin{enumerate}
  \item The normal fan $\Sigma_P$ of $P$ is complete.
        Whence every vector ${\bf v}$ can be written as a positive combination
        of the generators of a cone $\sigma \in \Sigma_P$.
  \item $\left( - \sum_{i=1}^n c_i {\bf v}^{j_i} \right)_k = \log |x_k|$,
        so ${\widetilde x}_k = (0.1)^{\log|x_k|} x_k = \frac{x_k}{|x_k|}$,
        whence $|{\widetilde x}_k| = 1$.

        $0 < z_{j_i} \leq 1 \Rightarrow 0 < {\widetilde z}_{j_i}
           = (0.1)^{c_i} z_{j_i} \leq 1$, as $c_i \geq 0$.
        Thus $0 < {\widetilde z}_j \leq 1$.~\qed
\end{enumerate}

\medskip

\noindent Note that any base $b$ can be used in the logarithm,
it suffices to replace 0.1 in Proposition~\ref{propscaling} by~$b^{-1}$.

\begin{example} (toric scaling analogue to projective space) {\rm
Suppose the Newton polytope is a simplex.  The homogeneous coordinates are
\begin{equation}
  \begin{tabular}{ccccccccccc}
      ( & $\lambda^1 x_1$    & , & $\lambda^0 x_2$ & ,
        & $\lambda^{-1} z_1$ & , & $z_2$ & , & $z_3$ & ) \\
      ( & $\lambda^0 x_1$    & , & $\lambda^1 x_2$ & ,
        & $z_1$ & , & $\lambda^{-1} z_2$ & , & $z_3$ & ) \\
      ( & $\lambda^{-1} x_1$ & , & $\lambda^{-1} x_2$ & ,
        & $z_1$ & , & $z_2$ & , & $\lambda^{-1} z_3$ & ) \\
  \end{tabular}
\end{equation}
for the respective facet normals ${\bf v}^1 = (1,0)$, ${\bf v}^2 = (0,1)$,
and ${\bf v}^3 = (-1,-1)$.
Since this is a {\em toric} analogue, we are concerned about solutions
with zero components and have more $z$-variables than in the case of
ordinary homogenization.  We consider two situations:
\begin{enumerate}
  \item $(x_1,x_2) = (2,2)$,
        ${\bf v} = (-\log(2),-\log(2))$,
        $J=(3,3)$, ${\bf c} = (\log(2),0)$
        \newline $\Rightarrow$ $z_3 = (0.1)^{\log(2)} z_3 = \frac{1}{2} z_3$.
  \item $(x_1,x_2) = (2,1)$,
        ${\bf v} = (-\log(2),-\log(1))$,
        $J=(2,3)$, ${\bf c} = (\log(2),\log(2))$
        \newline $\Rightarrow$ $z_2 = (0.1)^{\log(2)} z_2 = \frac{1}{2} z_2$
                           and $z_3 = (0.1)^{\log(2)} z_3 = \frac{1}{2} z_3$.
\end{enumerate}
In the first case it suffices to use one of the homogeneous coordinates,
while in the second case a combination is constructed.~\qex }
\end{example}

\begin{example} (Example~{\rm \ref{runningex}} continued.) {\rm
Here we also consider two situations:
\begin{enumerate}
  \item $(x_1,x_2) = (\frac{1}{2},2)$,
        ${\bf v} = (-\log(\frac{1}{2}),-\log(2))$,
        $J=(4,4)$, ${\bf c} = (\log(2),0)$
        \newline $\Rightarrow$ $z_4 = (0.1)^{\log(2)} z_4 = \frac{1}{2} z_4$.
  \item $(x_1,x_2) = (2,2)$,
        ${\bf v} = (-\log(2),-\log(1))$,
        $J=(2,3)$, ${\bf c} = (\log(2),\log(2))$
        \newline $\Rightarrow$ $z_3 = (0.1)^{2\log(2)} z_3 = \frac{1}{4} z_3$
                           and $z_4 = (0.1)^{3\log(2)} z_4 = \frac{1}{8} z_4$.
\end{enumerate}           
The scaling relation is constructed by decomposing $\bf v$ in $\Sigma_P$,
see Figure~\ref{figlatt}.~\qex}
\end{example}

In the general case when the system has an $r$-tuple of different
supports~$\cal A$, we decompose ${\bf v} = - \log|{\bf x}|$ in every
component of
$\Sigma_{\cal P} = (\Sigma_{P_1},\Sigma_{P_2}, \ldots, \Sigma_{P_r})$ to
compute new values of the $z$-variables which are distinct for every 
different support.

\smallskip

Propositions~\ref{propfacesols} and~\ref{propscaling} allow us to give a
concise proof of Bernshte\v{\i}n's second theorem formulated below as
Theorem~\ref{theoberntwo} using a homotopy argument.

\begin{definition}  {\rm
We say that the homotopy $H({\bf x},t) = {\bf 0}$ is {\em in generic position} 
with respect to the supports $\cal A$ that span the Newton polytopes $\cal P$,
if $\forall t \in [0,1)$: $H({\bf x},t) = {\bf 0}$ has supports~$\cal A$
and has exactly~$V_n({\cal P})$ regular solutions in~$(\cc^*)^n$. }
\end{definition}
Given a system~$F({\bf x}) = {\bf 0}$, the canonical example of
a homotopy in generic position is
\begin{equation}
  H({\bf x},t) = (1-t) F^{(0)}({\bf x}) + t F({\bf x}) = {\bf 0},
  \quad t \in [0,1],
\end{equation}
where the system~$F^{(0)}({\bf x}) = {\bf 0}$ has random complex
coefficients and the same supports as~$F$.

\begin{theorem} \label{theoberntwo}
Consider~$F({\bf x}) = {\bf 0}$, with supports~$\cal A$,
Newton polytopes~$\cal P$ and~$V_n({\cal P}) \not= 0$.
\begin{displaymath}
 \# \{ \ m_{\bf x} \ | \ F({\bf x}) = {\bf 0}, {\bf x}~{\rm isolated} \ \}
 < V_n({\cal P})
 ~ \Leftrightarrow ~
  \exists {\bf v} \in \rr^n \setminus \{ {\bf 0} \},
  \exists {\bf x} \in (\cc^*)^n:
  \partial_{\bf v} F({\bf x}) = {\bf 0}
\end{displaymath}
where $m_{\bf x}$ denotes the multiplicity of the
isolated solution $\bf x$ of $F({\bf x}) = {\bf 0}$.
\end{theorem}
{\em Proof}.  Construct any homotopy $H({\bf x},t) = {\bf 0}$ that is in
generic position with respect to~$\cal P$ and
for which $H({\bf x},t=1) = F({\bf x})$.  Consider the homogenized 
homotopy $\phi_{\cal A}(H)({\bf x},{\bf z},t) = {\bf 0}$.
\begin{itemize}
  \item[$\Rightarrow$] If $F({\bf x}) = {\bf 0}$ has fewer than~$V_n({\cal P})$
isolated solutions in~$(\cc^*)^n$, then, as $t \rightarrow 1$, either some
solution paths must diverge outside~$(\cc^*)^n$ or some solutions paths must
go to a component of solutions, because the homotopy is in generic position.
In the first case,
by Proposition~\ref{propscaling}, the only way to diverge is that some
$z_i \rightarrow 0$.  Then, according to Proposition~\ref{propfacesols}, 
we have a solution in $(\cc^*)^n$ of some face system.
In the second case, the solution component allows at least one free variable
to choose in moving on that component.  So we can move a free variable
outside $(\cc^*)^n$ for which some $z_i$ is zero.
Then we can again apply Proposition~\ref{propfacesols}.
  \item[$\Leftarrow$] If there is a face system that has a solution 
${\bf x} \in (\cc^*)^n$, then in the corresponding tuple of $z$-variables
for which $\phi_{\cal A}(H)({\bf x},{\bf z},1) = {\bf 0}$, 
$\exists i$: $z_i = 0$.
Since the homotopy is in generic position, there is at least one solution
path that ended at~$({\bf x},{\bf z})$.  By Proposition~\ref{propfacesols},
the transition to affine coordinates shows that the solution path
that ended at this couple~$({\bf x},{\bf z})$ diverged to a solution that
has components at infinity and/or zero components.
So $F({\bf x}) = {\bf 0}$ must have fewer than~$V_n({\cal P})$ isolated
solutions in~$(\cc^*)^n$.~\qed
\end{itemize}

\section{The Static Scaling in Projective Newton}

Let us recall Newton's method in projective space (see~\cite{li97}).
We assume we have a system $H = (h_1,h_2,\ldots,h_n)$
of homogeneous polynomials:
$h_i(\lambda {\bf z}) = \lambda^{\deg(h_i)} h_i({\bf z})$ for
$i= 1,2,\ldots,n$ with
$\lambda {\bf z} = (\lambda z_0, \lambda z_1, \ldots , \lambda z_n)$
the usual homogeneous coordinates.
For a solution at infinity: $z_0 = 0$.

\begin{algorithm} \label{algprojnewt}
(Newton's method in projective space) {\rm
%
\newline \begin{tabular}{lcr}
  Input: $(H,{\bf z})$. & \ \ \ \ 
				   & {\em (Homogeneous system, approximate root)} \\
  Output: $\bf z$.  & & {\em refined approximation} \\
\\
  ${\displaystyle i : |z_i| = \max_{j=1}^n |z_j|}$;
											   & & {\em maximal component} \\
  ${\bf z} := \left(
         \underbrace{\frac{z_0}{z_i}, \ldots, \frac{z_{i-1}}{z_i}} , 1 ,
         \underbrace{\frac{z_{i+1}}{z_i} , \ldots , \frac{z_n}{z_i}} \right)$;
                                        & & {\em scale the approximation} \\
  ${\bf x} := \left(
         ~ \overbrace{z_1,\ldots,z_{i-1}} ~ ,
         ~ \overbrace{z_{i+1},\ldots,z_n} ~
              \right)$;                  & & {\em switch to an affine chart} \\
  $F := \left. H \right|_{z_i = 1}$;   & & {\em select evaluation map} \\
  ${\bf x} := {\bf x} - D_{\bf x}^{-1} F({\bf x}) F({\bf x})$;
                                      & & {\em affine Newton method} \\
  ${\bf z} := (x_1,\ldots,x_i,1,x_{i+1},\ldots,x_n)$.
                                      & & {\em refined approximation} 
\end{tabular} }
\end{algorithm}
In the context of path following, $H$ is a parametrized family of systems.

\begin{example} (projective Newton on a sparse system) {\rm
Consider the application of projective Newton on a sparse system
of polynomials.  Suppose the supports of the polynomials are all equal.
We see the Newton polytope displayed in Figure~\ref{fignewt}.
In the homogenization the simplex is implicitly taken as Newton polytope
and to every facet we associate a new variable~$z_i$.
Since in projective space one does not distinguish solutions with
zero components, the original monomials are not copied.

\begin{figure}[hbt]
\begin{center}
\begin{picture}(420,100)(0,0)

\put(0,0){
\begin{picture}(120,100)(0,0)

\put(30,40){\circle*{3}}
\put(60,10){\circle*{3}}
\put(90,10){\circle*{3}}
\put(60,70){\circle*{3}}
\put(27,70){\line(1,0){6}}

\thicklines
\put(30,40){\line(1,-1){30}} 
\put(30,40){\line(1,1){30}} 
\put(60,10){\line(1,0){30}} 
\put(60,70){\line(1,-2){30}}

\setplotarea x from 0 to 60, y from 0 to 60
\setdashes <4pt>
\setlinear
\plot  30  10   30 100 /
\plot  30  10   120 10 /
\plot  30 100   120 10 /

\put(25,4){${}_0$}
\put(58,2){${}_1$}
\put(88,2){${}_2$}
\put(118,2){${}_3$}
\put(22,39){${}_1$}
\put(22,69){${}_2$}
\put(22,99){${}_3$}

\put(80,60){$z_0$}
\put(70,-3){$z_2$}
\put(13,55){$z_1$}

\put(50,35){$P$}

\end{picture}
}

\put(200,60){
$f({\bf x}) = c_{10} x_1 + c_{01} x_2 + c_{20} x_1^2 + c_{12} x_1 x_2^2$}

\put(200,30){
${\widetilde f}({\bf z}) 
= c_{10} z_0^2 z_1 + c_{01} z_0^2 z_2 + c_{20} z_0 z_1^2 + c_{12} z_1 z_2^2$}

\end{picture}

\caption{On the left we see how the polytope $P$ of $f$ is approximated
by a simplex which corresponds to the usual projective homogenization
${\widetilde f}$ given on the right.}
\label{fignewt}
\end{center}
\end{figure}

In case of divergence, $z_1$ and/or $z_2$ go to infinity.
The scaling step in Algorithm~\ref{algprojnewt} selects $i=1$ or $i=2$
and $z_0$ strives to zero.
In that case we see in Figure~\ref{fignewt} that the monomial
$z_1 z_2^2$ dominates.
However, a diverging path cannot go to a monomial and goes to one
of the other face systems.
The multi-homogeneous version of Newton's method will make the same
monomial dominant.~\qex }
\end{example}

\begin{example} (projective Newton leads to clustering) {\rm
Consider the following simple system:
\begin{equation}
   \left\{
      \begin{array}{r}
        x_1 - 10^8 = 0 \\
        x_2^2 - 1 = 0
      \end{array}
   \right.
\end{equation}
In affine space, the solutions $(x_1,x_2) = (10^8,\pm 1)$ are well-separated.
In projective space, they are represented as $(z_0,z_1,z_2) = (1,10^8,\pm 1)$
and scaled to $(z_0,z_1,z_2) = (10^{-8},1,\pm 10^{-8})$, which is a pair of
numerically clustered solutions.
Although this example is cooked up for this special occasion,
this phenomenon is observed as typical in many applications.
As the sum of the two Newton polytopes of the two polynomials is a rectangle,
the toric homogenization corresponds to the bi-homogenization, which
allows to scale the variables independently:
$(x_1,x_2) = (10^8,\pm 1)
\mapsto (z_{10},z_{11},z_{20},z_{21}) = (10^{-8},1,1,\pm 1)$,
avoiding the clustering of solutions.~\qex }
\end{example}

The most important difference with multi-projective Newton is the 
dynamic use of the scaling procedure of Proposition~\ref{propscaling}
which guarantees distinct paths, as formulated below in 
Theorem~\ref{theoseparate}.

\section{Toric Newton scales in a Dynamic Way}

The set up of Toric Newton enables an efficient numerical scaling
while path following.

\begin{algorithm} \label{algtorsymnewt} (Toric Newton -- symbolic set-up)
{\rm
%
\newline \begin{tabular}{lcr}
  Input: $(A,\triangle)$. & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
                          & {\em (support , triangulation)} \\
  Output: $H$;  & & {\em homogeneous coordinates matrix} \\
~~~~~~~~~~~ $\Sigma_P$. & & {\em simplicial normal fan} \\
\\
  1. $V$ := Facet\_Normals$(A,\triangle)$; & & {\em extract facet normals} \\
  2. $H$ := Homogenize$(A,V)$;        & & {\em relative lattice distances} \\
  3. $\Sigma_P$ := Cone\_Subdivision$(H,A,V)$. 
                               & & {\em refine normal fan}
\end{tabular} }
\end{algorithm}
Next we analyze the cost of this algorithm.
The dimension $n$ and the cardinality $\#A$ measure the input size.
The complexity of the output is bounded from below by the number of
cells $\# \triangle$ in any triangulation $\triangle$ of~$A$.
Recall that we are interested in {\em all} isolated solutions and to
count the number of roots we compute volumes.  So $\# \triangle$ is a
lower bound for the number of isolated solutions.

\begin{proposition} (Cost of Step 1 in Algorithm~\ref{algtorsymnewt}.)
Consider a triangulation $\triangle$ of the polytope $P$ spanned by~$A$.
Then all facets of~$P$ can be enumerated in
$O(\# \triangle n ( n^3 + \#A ))$ operations.
\end{proposition}
{\em Proof.}  The candidate facets are the facets of the simplices in
the triangulation.  Every simplex has $n+1$ facets.  The computation of
each candidate requires $O(n^3)$ arithmetical operations.
To test whether this candidate normal constitutes a facet normal we
have to compute inner products with the points in $A$ and look whether
the inner product is minimal for points on the facet of the simplex.
This gives the bound $O(\# \triangle n ( n^3 + \#A ))$.~\qed

\medskip

The number of arithmetic operations of Step 2 of Algorithm~\ref{algtorsymnewt}
equals~$O(n \# V)$.
Step 3 is only needed when the polytope is not simple, 
i.e.: when some vertices are contained in more than $n$ facets,
see the top of the polytope in Figure~\ref{figconeref}.
In Step~3 we make the fan simplicial subdividing the normal cones.
One method to subdivide normal cones is to cut corners
(see~\cite[page~190]{gkz94} and~\cite[pages~54--55]{zie95})
at vertices that have more than $n$ facets, as illustrated in 
Figure~\ref{figconeref}.
We refer to~\cite[Section~2.6]{ful93} for more on refinements of fans 
in connection with toric varieties.

\smallskip

For a normal cone $\sigma \in \Sigma_P$ at the vertex ${\bf a}$ we take the 
average~${\overline {\bf v}}$ of the facet normals that generate~$\sigma$.
Determine $c$ such that the hyperplane
$\langle {\bf x} , {\overline {\bf v}} \rangle = c$ separates ${\bf a}$
from the rest of the vertices that span the polytope.
The normal cones $\tau$ to the intersection points of the hyperplane 
with the edges of the polytope refine the normal cone~$\sigma$.
The recursive structure of the normal fan $\Sigma_P$ is illustrated at
the right of Figure~\ref{figconeref}.  In the decomposition of a vector
${\bf v}$ in a subdivision of $\sigma$, we replace the average
vector ${\overline {\bf v}}$ in terms of the other generators of $\sigma$,
so that we obtain a decomposition that uses only facet normals.

\begin{figure}[hbt]
\begin{center}
\begin{picture}(400,120)(0,0)

\put(-30,0){
\begin{picture}(120,120)(0,0)

\put(30,60){\vector(0,1){50}}    \put(18,110){$a_3$}
\put(30,60){\vector(1,0){100}}   \put(130,55){$a_2$}
\put(30,60){\vector(1,-4){15}}   \put(30,0){$a_1$}

\put(30,60){\circle*{3}}   \put(15,63){${}_{000}$}
\put(90,60){\circle*{3}}   \put(92,63){${}_{020}$}
\put(40,20){\circle*{3}}   \put(25,18){${}_{200}$}
\put(100,20){\circle*{3}}  \put(102,18){${}_{220}$}
% \put(65,40){\circle*{3}}        % 110
\put(65,95){\circle*{3}}   \put(59,100){${}_{111}$}

\put(45,85){{\scriptsize ${\bf b}^1$}}
\put(47,70){{\scriptsize ${\bf b}^2$}}
\put(75,85){{\scriptsize ${\bf b}^3$}}
\put(65,65){{\scriptsize ${\bf b}^4$}}

\thicklines
\put(30,60){\line(1,-4){10}}  
\put(90,60){\line(1,-4){10}}  
\put(30,60){\line(1,0){60}}  
\put(40,20){\line(1,0){60}}  

\setplotarea x from 0 to 120, y from 0 to 120
\setlinear
\plot   30 60  65 95 /
\plot   90 60  65 95 /
\plot   40 20  65 95 /
\plot  100 20  65 95 /

\put(55,85){\circle*{3}}
\put(58,73){\circle*{3}}
\put(75,73){\circle*{3}}
\put(72,85){\circle*{3}}

\plot  55 85  58 73 /
\plot  55 85  72 85 /
\plot  58 73  75 73 /
\plot  72 85  75 73 /

\end{picture}
}

\put(160,0){
\begin{picture}(200,120)(0,0)

\put(0,90){{\scriptsize $A = \left[ \begin{array}{c}
                           1~1~1 \\
                           0~0~0 \\
                           2~0~0 \\
                           0~2~0 \\
                           2~2~0 \\
                         \end{array}
                 \right]$}}

\put(70,90){{\scriptsize $V = \left[ \begin{array}{r}
                            0~\!~~~~0~\!~~~~1 \\
                            0~\!~~~~1~-1 \\
                           -1~\!~~~~0~-1 \\
                            1~\!~~~~0~-1 \\
                            0~-1~-1 \\
                         \end{array}
                 \right]$}}

\put(155,110){${\bf v}^1$}
\put(155,100){${\bf v}^2$}
\put(155,90){${\bf v}^3$}
\put(155,80){${\bf v}^4$}
\put(155,70){${\bf v}^5$}

\put(190,90){{\scriptsize $H = \left[ \begin{array}{c}
                           1~0~0~0~0 \\
                           0~0~2~0~2 \\
                           0~0~0~2~2 \\
                           0~2~2~0~0 \\
                           0~2~0~2~0 \\
                         \end{array}
                 \right]$}}

\put(0,30){{\scriptsize $\Sigma_P = \left[ \begin{array}{c}
                           1~1~1 \\
                           0~0~0 \\
                           2~0~0 \\
                           0~2~0 \\
                           2~2~0 \\
                         \end{array}
                 \right]$}}

\put(57,50){\vector(1,0){100}}
\put(57,41){\vector(1,0){20}}
\put(57,32){\vector(1,0){20}}
\put(57,23){\vector(1,0){20}}
\put(57,14){\vector(1,0){20}}

\put(80,10){$({\bf v}^1,{\bf v}^3,{\bf v}^5)$}
\put(80,19){$({\bf v}^1,{\bf v}^4,{\bf v}^5)$}
\put(80,28){$({\bf v}^1,{\bf v}^2,{\bf v}^3)$}
\put(80,37){$({\bf v}^1,{\bf v}^2,{\bf v}^4)$}

\put(160,30){{\scriptsize $\left[ \begin{array}{c}
                                    {\bf b}^1 \\
                                    {\bf b}^2 \\
                                    {\bf b}^3 \\
                                    {\bf b}^4 \\
                                 \end{array}
                           \right]$}}

\put(185,45){\vector(1,0){20}}
\put(185,36){\vector(1,0){20}}
\put(185,27){\vector(1,0){20}}
\put(185,18){\vector(1,0){20}}

\put(210,41){$({\overline {\bf v}},{\bf v}^2,{\bf v}^4)$}
\put(210,32){$({\overline {\bf v}},{\bf v}^2,{\bf v}^3)$}
\put(210,23){$({\overline {\bf v}},{\bf v}^4,{\bf v}^5)$}
\put(210,14){$({\overline {\bf v}},{\bf v}^3,{\bf v}^5)$}

\end{picture}
}

\end{picture}

\caption{Refinement of a normal cone of a nonsimple polytope by cutting the
corner at~$(1,1,1)$.  The points ${\bf b}^i$ lie on the intersection of the
hyperplane cut with the edges that contain $(1,1,1)$. }
\label{figconeref}
\end{center}
\end{figure}

In three dimensions we use only one separating hyperplane,
because every edge is enclosed in exactly two facets.
This is no longer the case in higher dimensions, so we proceed
recursively as the points in the intersection of the polytope with
the separating hyperplane belong to a space of lower dimension.
By Proposition~2.4 from~\cite[page~55]{zie95}, edges are in one-to-one
correspondence with vertices of the polytope cut out by the separating
hyperplane.  Recursive application proves the following.

\smallskip

\begin{proposition} \label{propcostconeref}
(Cost of Step 3 in Algorithm~\ref{algtorsymnewt}.) 
To triangulate an $n$-dimensional cone at most $n-2$ separating hyperplanes
are needed.
\end{proposition}
Although the number of cones grows exponentially, one new
refining cone requires only elementary operations to compute the 
intersection points of the edges with the separating hyperplane.
Therefore we can perform a lazy refinement of $\Sigma_P$, decomposing
only those normal cones when needed in the scaling procedure.

\medskip

Just as projective Newton method, we switch to an appropriate
affine chart of dimension~$n$.

\begin{algorithm} \label{algtornumnewt}
(Toric Newton -- number crunching) {\rm
\newline \begin{tabular}{lcr}
  Input: $\phi = \phi_A(F)$; & \ \ \ \ \ \ \ \ \ \ \
                              & {\em system in homogeneous coordinates} \\
~~~~~~~~~ $({\bf x},{\bf z})$;
  & & {\em approximate root with $|x_i| = 1$ and $|z_j| \leq 1$} \\
~~~~~~~~~ $\Sigma_P$. & & {\em simplicial normal fan} \\
  Output: $({\bf x},{\bf z})$.
  & & {\em refined approximation with $|x_i| = 1$ and $|z_j| \leq 1$} \\
\end{tabular}

\begin{tabular}{lcr}
  $F := \left. \phi \right|_{Fix~{\bf z}}$; & & {\em evaluation map} \\
  ${\bf x} := {\bf x} - D_{\bf x}^{-1} F({\bf x}) F({\bf x})$;
           & \ \ \ \ \ \ \ \ \ \ \ \ & {\em affine Newton method} \\  
  ${\bf v} := -\log|{\bf x}|$;  & & {\em direction of correction} \\
  $J = (j_1,j_2,\ldots,j_n)$
  & &  {\em find $\sigma \in \Sigma_P : {\bf v} \in \sigma$} \\
  ~~~ := Localize\_Position$(\Sigma_P,{\bf v})$;
  & & {\em $\sigma = 
       {\rm span}({\bf v}^{j_1},{\bf v}^{j_2},\ldots,{\bf v}^{j_n})$} \\
  ${\displaystyle {\bf c} = (c_1,c_2,\ldots,c_n) :
      {\bf v} = \sum_{i=1}^n c_i {\bf v}^{j_i}}$;
  & & {\em decompose ${\bf v}$} \\
  for all $i \in \{1,2,\ldots,n \}$ do & & \\
  ~~ $x_i := \frac{x_i}{|x_i|}$; $z_{j_i} := (0.1)^{c_i} z_{j_i}$;
  & & {\em scale solution and coefficients} \\
  end for. & & 
\end{tabular} }
\end{algorithm}

In analyzing the cost of localizing the position of the direction vector
${\bf v}$ we assume that the fan $\Sigma_P$ has a recursive simplicial
structure, as illustrated in Figure~\ref{figconeref}.

\begin{proposition} (Cost of {\rm Localize\_Position} in
Algorithm~\ref{algtornumnewt}.)
For ${\bf v} \in \rr^n$ and $\Sigma_P$, $P = {\rm conv}(A)$,
finding $\sigma \in \Sigma_P$ so that ${\bf v} \in \sigma$ takes 
at least $\#A$ and at most $\#A (n-1)$ inner products.
\end{proposition}
{\em Proof.}  If $P$ is simple, then $\sigma$ is the normal cone
at the point ${\bf a} \in A$ for
which $\langle {\bf a} , {\bf v} \rangle
= \min_{{\bf b} \in A} \langle {\bf b} , {\bf v} \rangle$.
Computing this point requires $\#A$ inner products.
If $P$ is not simple, there are at most $n-2$ refinement levels,
according to Proposition~\ref{propcostconeref}.
At every level, there are at most $\#A$ points to consider.~\qed

\medskip

For $r$-tuples of supports~$\cal A$, we construct
$\Sigma_{\cal P} = (\Sigma_{P_1},\Sigma_{P_2}, \ldots, \Sigma_{P_r})$.
The cost of the geometric constructions is scaled by a factor~$r$.
  From our experiences with mixed-volume computations~\cite{vgc96}
we find linear upscaling in the cost rather remarkable.

\section{Numerical Convergence Properties}

In the absence of singularities, Newton's method converges quadratically,
doubling the precision of the approximate root in each step.
With a homotopy in generic position, we will not hit a solution
component, except perhaps only at the very end of the paths.
What may happen however, is that paths become clustered, hereby
not only limiting the convergence but possibly also terminating
the numerical path tracking with a failure to reach the end.

\smallskip

Affine path following methods deal well with paths diverging
outside $(\cc^*)^n$, for those paths exhibit a steady direction leading
to the face systems that cause the deficiency.  To obtain a certificate
of divergence, polyhedral end games~\cite{hv98} can effectively trace
the paths in affine space, because to compute the direction of divergence
one does not need to follow the paths up to infinity.
Reaching badly-scaled though probably still meaningful solutions by
affine path following methods is often impossible due to the finiteness
of numbers a computer can represent.  Projective Newton brings infinity
to our working space, but roots that are clearly separated in affine space
can cluster by a projective transformation.

\smallskip

The dynamic of the scaling procedure of Proposition~\ref{propscaling} provides
the proof of the numerical convergence properties of Toric Newton.

\smallskip

\begin{theorem} \label{theoseparate}
The solution paths defined by a homotopy in generic position compactified
torically remain separated, except possibly at the end when converging
to multiple solutions.
\end{theorem}
{\em Proof}.  At the start of the paths the solutions $({\bf x},{\bf z})$
are all mutually distinct tuples.  We will show that they stay separated,
as the continuation parameter $t$ progresses.
By Proposition~\ref{propscaling}, the components of ${\bf x}$ remain on
the unit circle and the components of ${\bf z}$ may possibly go from 1 to~0.

\smallskip

Either a path stays inside $(\cc^*)^n$ or diverges, as $t \rightarrow 1$.
Only in the latter case, according to Theorem~\ref{theoberntwo}, the
solution paths have the same direction.
So if the path stays inside $(\cc^*)^n$, its direction ${\bf v}$ 
in Proposition~\ref{propscaling} is different from the other paths
yielding distinct values for ${\bf z}$.
Thus all paths that stay inside $(\cc^*)^n$ remain separated, unless they
converge to a multiple solution.

\smallskip

Still according to Theorem~\ref{theoberntwo},
only paths that diverge to solutions of the 
same $\partial_{\bf v} F({\bf x}) = {\bf 0}$ can have
the same direction~${\bf v}$.  However, this direction becomes only apparent
at the end of the paths, when the monomials of
$\partial_{\bf v} F({\bf x}) = {\bf 0}$ are dominating the rest of the system.
In that case, even when two paths would converge to the same ${\bf z}$-vector,
they would converge to the same ${\bf x}$-vector only if that vector
would be a multiple solution of~$\partial_{\bf v} F({\bf x}) = {\bf 0}$.~\qed

\medskip

The following is an immediate reformulation of Theorem~\ref{theoseparate}.

\begin{corollary}
For a homotopy in generic position, Toric Newton method converges
quadratically in correcting the solutions along the paths, except
when approaching singularities at the end.
Regular solutions are approximated by a quadratically
converging Toric Newton method.
\end{corollary}

We close this section with an important numerical observation.
One may object that as $z$-variables can go to zero we could have
to rely on very tiny values to distinguish the solutions.
However, we must not look at the actual values of the $z$-variables,
but at their exponents.  These exponents have been obtained in a
numerically stable way, by adding consecutive coefficients in the
decomposition of the directions of the paths, see the scaling procedure of
Proposition~\ref{propscaling}.
To deal with extremal values of solutions, we usually work with a
logarithmic scale, as the magnitude of the solutions is most relevant.

\section{Computational Experiences}

Here we report on the geometric computations showing
the practical feasibility of the method.

\subsection{Implementation Issues}

The experiments are conducted with the polynomial homotopy continuation
environment~PHC~\cite{ver97} created, developed and maintained by the author.
On the occasion of reported calculations with SAGBI homotopies~\cite{ver98},
a new version of the package was recently released.

\smallskip

The construction of the triangulation is performed
with the dynamic lifting algorithm~\cite{vgc96}.
This algorithm is efficient for sparse inputs.
To overcome the space complexity, one should refine the fan 
in a lazy way subdividing
only those parts of the cones that are strictly needed while path following.
The implementation of Toric Newton produces values 
$(\theta,{\bf z}) \in [0,2 \pi)^n \times [0,1]^N$, where
$x_i = e^{i \theta_i}$.  To minimize numerical round off, the affine
Newton's method must be executed in~$[0,2 \pi)^n \times [0,1]^N$.
This stage of implementation in~PHC is still under construction.

\smallskip

Multi-precision arithmetic (as available in the latest release of~PHC)
is recommended for transforming from toric to affine
coordinates at the end of the paths and for evaluation of the residuals.
Only and only at that stage, we fine multi-precision arithmetic useful.
Alternatively and numerically, one could treat the $z_i$'s one after the
other and apply Newton's method intermediately to correct the rounding errors.
For well-conditioned affine roots, this method works well.

\subsection{Geometric Samples}

Here we give an impression on the complexity of the geometry involved
on a 15-point configuration randomly generated for several dimensions.
Table~\ref{tabgeosamples} summarizes the characteristics and timings.
The computations were done on a
166 MHz Pentium II processor with 64 Mb internal memory running Linux.

{\small
\begin{table}[ht]
\begin{center}
\begin{tabular}{|c||r|r|r||r|r|r|r||r|r|} \hline
  & \multicolumn{7}{c||}{15 points in $\{ 0,1,2,3 \}^n$ for $n=3,4,5,6,7$}
  & \multicolumn{2}{c|}{cpu min, sec, millisec} \\
 \cline{2-10} 
n & $\# \triangle$ & vol & $\# V$ & max $\omega$ &
    ${\scriptsize \deg(\phi_A,z_i)}$ 
  & $\# \sigma$ & $\# \Sigma_P$ & construct & decompose \\ \hline \hline
3 &  16 &     64 &  17 &         57 &    15 &  6 &     63
  &      20~ms & 380~ms \\
4 &  21 &    172 &  25 &      1,227 &    36 & 15 &    341
  &      90~ms & 500~ms \\
5 &  57 &  1,081 &  59 &      5,904 &    70 & 27 & 17,299
  &  5~s~40~ms & 850~ms \\
6 & 110 &  4,159 & 140 &  3,619,556 &   660 & 78 & 55,725
  & 20~s~30~ms & 1~s~450~ms \\ 
7 & 157 & 29,041 & 248 & 12,567,604 & 1,364 & 135 & 664,135
  & 5~m~31~s~470~ms & 3~s~~30~ms \\ \hline 
\end{tabular}
\caption{Characteristics and timings for random samples.
The first part lists the number of cells~$\# \triangle$, the volume of the
polytope and the number of facet normals~$\# V$.  In the second part we 
compare the maximal height (max~$\omega$) to construct~$\triangle$ to
the highest degree of a $z$-variable in the homogenization.
The sixth column displays the maximum number of generating facet 
normals~$\# \sigma$ of a cone and the total number of cones in the
refinement of the fan~$\Sigma_P$.  The last part shows timings for
constructing this refinement and for decomposing 1000 random vectors.}
\label{tabgeosamples}
\end{center}
\end{table}
}

  From Table~\ref{tabgeosamples} we observe the growth in complexity as
the dimension increases.
The maximal lifting value~$\omega$ is generated by the dynamic lifting
algorithm~\cite{vgc96}.
Although the integer lifting values generated by the implementation
of the dynamic lifting algorithm can be improved by balanced floating-point
liftings developed in~\cite{glvw97}, it is noteworthy that the degrees of
the $z$-variables are of a much lesser magnitude.

\smallskip

Newton polytopes are in general not simple, as experienced from
the~$\# \sigma$.
The exhaustive refinement of the normal fan is quite memory exhaustive.
For larger dimensions we must use a lazy decomposition of normal cones.
The fast timings in constructing and consulting this geometric data
structure are promising and encouraging.

\subsection{A Numerical Example}

The following system is inspired on an example in~\cite{kp96}.
\begin{displaymath}
  \left\{
    \begin{array}{c}
       x_1 - 10 = 0 \\
       x_2 - x_1^2 = 0 \\
       x_3 - x_2^2 = 0 \\
       x_4 - x_3^2 = 0 \\
    \end{array}
  \right.
  \quad \Leftrightarrow \quad
  \left\{
    \begin{array}{c}
       \gamma_{11} (x_1 - 10) + \gamma_{12} (x_2 - x_1^2)
       + \gamma_{13} (x_3 - x_2^2) + \gamma_{14} (x_4 - x_3^2) = 0 \\
       \gamma_{21} (x_1 - 10) + \gamma_{22} (x_2 - x_1^2)
       + \gamma_{23} (x_3 - x_2^2) + \gamma_{24} (x_4 - x_3^2) = 0 \\
       \gamma_{31} (x_1 - 10) + \gamma_{32} (x_2 - x_1^2)
       + \gamma_{33} (x_3 - x_2^2) + \gamma_{34} (x_4 - x_3^2) = 0 \\
       \gamma_{41} (x_1 - 10) + \gamma_{42} (x_2 - x_1^2)
       + \gamma_{43} (x_3 - x_2^2) + \gamma_{44} (x_4 - x_3^2) = 0 \\
    \end{array}
  \right.
\end{displaymath}
where the $\gamma$-coefficients are randomly generated complex numbers.
Both systems are equivalent.  The embedding at
the right is only meant to illustrate how four nice quadratic polynomial
equations can hide one root with large components.
This problem is not recognized by conventional numerical scaling methods
and is the kind of target application for Toric Newton.

\section{Closing Remarks}

To solve a wide class of polynomial systems in a reliable and efficient
manner by homotopy continuation, a variety of different methods has been
developed.  Coefficient and equation scaling, often used in combination
with centering, is a remedy to deal with systems that have extremal
coefficients.  Projective transformations compactify our working space.
Polyhedral end games used in affine path following methods certify
diverging paths.  Toric Newton generates an optimal affine chart to 
approximate the magnitude of extremal roots.

\smallskip

Much is left undone.  At a practical level, path trackers have to be
developed and implemented that recognize the toric situation.
Instead of complex condition numbers there is a need for structured
condition numbers that measure the accuracy of the computations.

\bigskip

\noindent {\bf Acknowledgments}.  This work could not have been started 
without Birk Huber with whom I first learned in~\cite{hv98} to compute
mathematical certificates of divergence.
The practical experiences on path following in projective space shared
to me by Tangan Gao were very inspiring.
I thank T.Y. Li for a hearty conversation that led to the inclusion of
assumption that homotopies need to be in generic position.
I am grateful for the conversations with Alicia Dickenstein and 
Eduardo Cattani.  

{\small
\begin{thebibliography}{10}

\bibitem{ber75}
D.N. Bernshte{\v{\i}}n.
\newblock The number of roots of a system of equations.
\newblock {\em Functional Anal. Appl.}, 9(3):183--185, 1975.
\newblock Translated from {\em Funktsional. Anal. i Prilozhen.},
  9(3):1--4,1975.

\bibitem{bcss97}
L.~Blum, F.~Cucker, M.~Shub, and S.~Smale.
\newblock {\em Complexity and Real Computation}.
\newblock Springer--Verlag, New York, 1997.

\bibitem{cd97}
E.~Cattani and A.~Dickenstein.
\newblock A global view of residues in the torus.
\newblock {\em Journal of Pure and Applied Algebra}, 117--118:119--144, 1997.
\newblock Proceedings of MEGA'96.

\bibitem{cds96}
E.~Cattani, A.~Dickenstein, and B.~Sturmfels.
\newblock Computing multidimensional residues.
\newblock In L.~Gonzalez-Vega and T.~Recio, editors, {\em Algorithms in
  Algebraic Geometry and Applications}, volume 143 of {\em Progress in
  Mathematics}, pages 135--164. Birkh\"auser, Basel, 1996.
\newblock Proceedings of MEGA'94.

\bibitem{cls98}
D.~Cox, J.~Little, and D.~O'Shea.
\newblock {\em Using Algebraic Geometry}, volume 185 of {\em Graduate Texts in
  Mathematics}.
\newblock Springer--Verlag, New York, 1998.

\bibitem{cox95}
D.A. Cox.
\newblock The homogeneous coordinate ring of a toric variety.
\newblock {\em J. Algebraic Geometry}, 4(1):17--50, 1995.

\bibitem{dan78}
V.~Danilov.
\newblock The geometry of toric varieties.
\newblock {\em Russ. Math. Surveys}, 33:97--154, 1978.

\bibitem{ded97}
J.P. Dedieu.
\newblock Condition number analysis for sparse polynomial systems.
\newblock In F.~Cucker and M.~Shub, editors, {\em Foundations of Computational
  Mathematics. Selected Papers of a Conference, Held at IMPA in Rio de Janeiro,
  January 1997}, pages 75--101. Springer--Verlag, 1997.

\bibitem{ds97}
J.P. Dedieu and M.~Shub.
\newblock Multihomogeneous {N}ewton methods.
\newblock {IBM} {T}echnical {R}eport, IBM Research Division, 1997.

\bibitem{ds98}
J.P. Dedieu and M.~Shub.
\newblock Newton and predictor-corrector methods for over-determined systems of
  equations.
\newblock {IBM} {T}echnical {R}eport, IBM Research Division, 1998.

\bibitem{ful93}
W.~Fulton.
\newblock {\em Introduction to Toric Varieties}, volume 131 of {\em Annals of
  Mathematics Studies}.
\newblock Princeton University Press, Princeton, New Jersey, 1993.

\bibitem{glvw97}
T.~Gao, T.Y. Li, J~Verschelde, and M.~Wu.
\newblock Balancing the lifting values to improve the numerical stability of
  polyhedral homotopy continuation methods.
\newblock Manuscript, 1997.

\bibitem{gkz94}
I.M. {Gel'fand}, M.M. Kapranov, and A.V. Zelevinsky.
\newblock {\em {D}iscriminants, {R}esultants and {M}ultidimensional
  {D}eterminants}.
\newblock Birkh\"auser, Basel, 1994.

\bibitem{hs95}
B.~Huber and B.~Sturmfels.
\newblock A polyhedral method for solving sparse polynomial systems.
\newblock {\em Math. Comp.}, 64(212):1541--1555, 1995.

\bibitem{hv98}
B.~Huber and J.~Verschelde.
\newblock Polyhedral end games for polynomial continuation.
\newblock {\em Numerical Algorithms}, 18(1):91--108, 1998.

\bibitem{kp96}
T.~Krick and L.M. Pardo.
\newblock A computational method for diophantine approximation.
\newblock In L.~Gonzalez-Vega and T.~Recio, editors, {\em Algorithms in
  Algebraic Geometry and Applications}, volume 143 of {\em Progress in
  Mathematics}, pages 193--253. Birkh\"auser, Basel, 1996.
\newblock Proceedings of MEGA'94.

\bibitem{li97}
T.Y. Li.
\newblock Numerical solution of multivariate polynomial systems by homotopy
  continuation methods.
\newblock {\em Acta Numerica}, 6:399--436, 1997.

\bibitem{mm87}
K.~Meintjes and A.P. Morgan.
\newblock A methodology for solving chemical equilibrium systems.
\newblock {\em Appl. Math. Comput.}, 22(4):333--361, 1987.

\bibitem{mor87}
A.~Morgan.
\newblock {\em Solving polynomial systems using continuation for engineering
  and scientific problems}.
\newblock Prentice-Hall, Englewood Cliffs, N.J., 1987.

\bibitem{oda88}
T.~Oda.
\newblock {\em Convex {B}odies and {A}lgebraic {G}eometry: an {I}ntroduction to
  the {T}heory of {T}oric {V}arieties}.
\newblock Springer--Verlag, Berlin Heidelberg New York, 1988.

\bibitem{roj98}
J.M. Rojas.
\newblock Solving degenerate sparse polynomial systems faster.
\newblock Manuscript available at http://www.cityu.edu.hk/ma/staff/rojas/.

\bibitem{rs98}
M.F. Roy and A.~Szpirglas.
\newblock Bezoutiens et r{\'e}sidus, affines, projectifs et dans le tore.
\newblock Manuscript, 1998.

\bibitem{ver97}
J.~Verschelde.
\newblock {PHCPACK}: A general-purpose solver for polynomial systems by
  homotopy continuation.
\newblock Paper and software available at the author's web-pages.

\bibitem{ver98}
J.~Verschelde.
\newblock Numerical evidence for a conjecture in real algebraic geometry.
\newblock Preprint \#1998-064, MSRI, 1998.
\newblock Paper and software available at the author's web-pages.

\bibitem{vgc96}
J.~Verschelde, K.~Gatermann, and R.~Cools.
\newblock Mixed-volume computation by dynamic lifting applied to polynomial
  system solving.
\newblock {\em Discrete Comput. Geom.}, 16(1):69--112, 1996.

\bibitem{vvc94}
J.~Verschelde, P.~Verlinden, and R.~Cools.
\newblock Homotopies exploiting {N}ewton polytopes for solving sparse
  polynomial systems.
\newblock {\em SIAM J. Numer. Anal.}, 31(3):915--930, 1994.

\bibitem{wbm87}
L.T. Watson, S.C. Billups, and A.P. Morgan.
\newblock Algorithm 652: {HOMPACK}: a suite of codes for globally convergent
  homotopy algorithms.
\newblock {\em ACM Trans. Math. Softw.}, 13(3):281--310, 1987.

\bibitem{wsmmw97}
L.T. Watson, M.~Sosonkina, R.C. Melville, A.P. Morgan, and H.F. Walker.
\newblock {HOMPACK90}: A suite of {F}ortran 90 codes for globally convergent
  homotopy algorithms.
\newblock {\em ACM Trans. Math. Softw.}, 23(4):514--549, 1997.
\newblock Available at http://www.cs.vt.edu/{\~{}}ltw/.

\bibitem{zie95}
G.M. Ziegler.
\newblock {\em Lectures on {P}olytopes}, volume 152 of {\em Graduate Texts in
  Mathematics}.
\newblock Springer--Verlag, New York, 1995.

\end{thebibliography}
}

\end{document}
