%\documentstyle{article}
%Uncomment the previous line and comment out the following two for latex 2.09 compatibility
\documentclass{article}
\usepackage{latexsym}

\begin{document}
\title{Numerical Computation of Polynomial Roots Using MPSolve \\
        Version  2.2}
\author{Dario Andrea Bini and Giuseppe Fiorentino \\
 Dipartimento di Matematica, Universit\`a di Pisa \\
 Via Buonarroti 2, 56127 Pisa. \\
 (bini@dm.unipi.it, fiorent@di.unipi.it)
}
\date{January 2000}
\maketitle

\begin{abstract}
  We describe the implementation and the use of the package MPSolve
  (Multiprecision Polynomial Solver) for the approximation of the roots
  of a univariate polynomial $p(x)=\sum_{i=0}^n a_i x^i$.
  
  The package relies on an algorithm, based on simultaneous
  approximation techniques, that generates a sequence ${\cal A}_1
  \supset {\cal A}_2 \supset \ldots \supset{\cal A}_k$ of nested sets
  containing the root neighborhood of $p(x)$ (defined by the precision
  of the coefficients) and outputs
  Newton-isolated approximations of the roots.  The algorithm,
  particularly suited to deal with sparse polynomials or polynomials
  defined by straight line programs, adaptively adjusts the working
  precision to the conditioning of each root. In this way, only the
  leading digits of the input coefficients sufficient to recover the
  requested information on the roots are actually involved in the
  computation.  This feature makes our algorithm particularly suited
  to deal with polynomials arising from certain symbolic computations
  where the coefficients are typically integers with hundreds or
  thousands digits.
  
  The present release of the software may perform different computations
  such as counting, isolating or approximating (to any number of
  digits)  the roots belonging to a given subset of the complex plane.
  Automatic detection of multiplicities and/or of real/imaginary roots
  can be also performed. Polynomials with approximately known
  coefficients are also allowed.
\end{abstract}

\section{Introduction}
MPSolve is a software tool for the numerical computation of the roots
of a univariate polynomial $p(x)=\sum_{i=0}^n a_i x^i$.
The executable program of the package is called {\tt unisolve} 
(univariate polynomial solver).

In the present release the following features are available.

\begin{description}
\item[Input:] \ \\

  \begin{itemize}
  \item The input polynomial can be provided either by means of the
    degree $n$ and the coefficients $a_i$, $i=0,\ldots, n$ of the
    univariate polynomial $p(x)$, or as a C program for the evaluation
    of the values of $p(x)$, of $p'(x)$ and of an upper bound of
    $|p(x)-{\rm fl}(p(x))|$, where fl($a$) denotes the value of the
    expression $a$ computed in floating point arithmetic.
   
 \item The user can define both the input precision $d_{in}$ and the
   desired output precision $d_{out}$ expressed as decimal digits.
   The input precision defines the {\it polynomial neighborhood} of
   $p(x)$, i.e., the set of all the polynomials with coefficients
   having $d_{in}$ common digits with the corresponding coefficients
   of $p(x)$, and the {\it root neighborhood}, i.e., the set of the
   roots of all the polynomials in the polynomial neighborhood of
   $p(x)$.  For polynomials having coefficients with infinite
   precision (say, integer or rational) $d_{in}$ must be set to zero.
  \end{itemize}

\item[Output:] \ \\
   In order to provide the user as much information as possible, MPsolve
   outputs its approximations according to the following rules:
   \begin{itemize}
   \item For infinite precision input, i.e., for $d_{in}=0$, the
     program delivers a list of $n$ complex numbers represented with
     at most $d_{out}$ digits.  For each number of this list, the
     displayed digits coincide with the digits of the corresponding
     root of $p(x)$. The number of displayed digits reaches (or even
     may slightly exceed) the maximum value $d_{out}$ only for those
     sets of roots that have at least $d_{out}$ common digits.
     Otherwise it is displayed a number of digits sufficient to
     provide Newton-isolated approximations of the roots.  That is, the
     number of digits is large enough to separate each root from the
     others and to guarantee the quadratic convergence of Newton's
     iteration right from the start when applied to the output
     approximation. For instance, for the polynomial having the
     following 5 roots 1.11111111, 1.12222222, 1.1233333, 1/3, 1/3,
     and for the input $d_{in}=0$, $d_{out}=30$ the program would
     output the numbers (1.111, 0.0) (1.1222, 0.0) (1.12333, 0.0)
     (0.333333333333333333333333333333, 0.0) \\
     (0.333333333333333333333333333333, 0.0).

   \item
   For input with a finite precision ($d_{in}>0$) the program delivers
   a list of $n$ complex numbers having at most $d_{out}$ digits as in
   the case of infinite precision. For each number $z$ in this list
   there exists a polynomial $q(x)$ in the neighborhood of $p(x)$ 
   such that $q(z)=0$. Each root of $p(x)$ has its corresponding
   approximation in the list. Moreover, any polynomial $q(x)$ in the
 neighborhood of $p(x)$ has a root that coincides with $z$ in the 
leading displayed digits.

   \item 
   If the output format -Of is chosen, the program delivers the
   numerical approximations to the roots with no filtering of the
   correct digits, a guaranteed a posteriori error bound, a string of
   three characters describing the status of the approximation.  This
   string has the following meaning: the first character may take the
   following values: 
   \begin{itemize} 
     \item m:  multiple root    
     \item i:  isolated root 
     \item a:  approximated single root (relative error less than 
                $10^{-d_{out}}$)
     \item o:  approximated cluster of roots (relative error less than
                $10^{-d_{out}}$)
     \item c:  cluster of roots not yet approximated (relative error
                greater than $10^{-d_{out}}$) 
   \end{itemize} 
   The second character may take the following
   values 
   \begin{itemize} 
       \item R:  real root
       \item r:  non-real root
       \item I:  imaginary root
       \item i:  non-imaginary root
       \item w:  uncertain real/imaginary root
       \item z:  non-real and non-imaginary root
   \end{itemize} 
   The third character may take the following values
   \begin{itemize}
       \item i:  root in $\cal S$
       \item o:  root out of  $\cal S$
       \item u:  root  uncertain
   \end{itemize}
 
   
 \item In certain cases it may happen that the relative error bound of
   the computed approximations is greater than 1, i.e., there is no
   correct bit in the output. This situation typically occurs, when
   the imaginary (or real) part of a nonzero root is zero or is such
   that its ratio with the real (imaginary) part has modulus less than
   $10^{-d_{out}}$. In this case the program outputs 0.0e{\em xx},
   where {\em xx} is the exponent of the approximation of this number.
   Therefore, an output like (1.234, 0.0e-1000) means that the
   imaginary part, if nonzero, has a modulus less than $10^{-1000}$.
   This lack of information is typical in numerical computation when
   we have to deal with 0. Moreover, for polynomials having
   coefficients with infinite precision, it is easy to overcome this
   lack of information by using the reality or integrality of the
   coefficients. This is performed automatically with the option {\tt
     -Dx}, where `{\tt x}' may be {\tt r} (real detect) {\tt i}
   (imaginary detect) or {\tt b} (both real and imaginary detect). See
   the next section about the available options.

   The uncertainty about the zero output is present also in the case
   of polynomials $p(x)$ having approximate input coefficients where
   the root neighborhood of $p(x)$ intersects
   zero. In this case the program outputs 0.0 for those approximations
   belonging to the connected component of the root neighborhood that
   intersects zero. However, these (nonzero) approximations are
   reported integrally in the standard error if the switch {\tt -d }
   (debug) is set or the output format {\tt -Of} is chosen 
   (see the next section for more details).
  \end{itemize}

\item[Goals:] \ \\
  The program has been designed to deal with different goals.  In
  particular besides the goal ``Newton isolation'' described above
  (default option {\tt -Gi}), the program allows the further goals
  ``approximation'' (option {\tt -Ga }) and ``count roots'' (option {\tt
  -Gc}).  Moreover, by giving some command line options, 
   the user is enabled to restrict the search of the
  roots to specific subsets of the complex plane (option {\tt -S}{\em
  x}), say, unit disk, right half-plane, etc.), also the automatic
  detection of multiple roots (option {\tt -M}) and/or of real and
  imaginary roots (option {\tt -D}) is fully implemented.

\end{description}

The program, based on the technique of simultaneous approximation of
the roots, makes use of a variant of the algorithms
described in [1], [2]. Essentially, the program constructs a sequence
of sets ${\cal A}_1 \supset {\cal A}_2 \supset {\cal A}_3 \supset
\ldots \supset{\cal A}_k$, where each set ${\cal A}_i$ is the union of
$n$ disks, contains the root neighborhood of $p(x)$, and each disk
contains at least one root. Moreover, each connected component of
${\cal A}_i$ contains as many roots of any polynomial in the
neighborhood of $p(x)$ as the number of disks that make up the
component.  Each set ${\cal A}_i$ is obtained by performing the
computation with increasing levels of working precision $w_1$, $w_2,
\ldots$, starting with the standard machine precision ($w_1=53$ bits
for the IEEE standard) and doubling the precision at each stage,
i.e. $w_{i+1}=2w_i$, until either all the disks are Newton--isolated
or the output precision has been reached, or the working precision
becomes larger than the input precision, i.e., $w_i>4 n d_{in}\log_2
10$.

The main engine for shrinking the disks is Aberth's iteration.  In
order to accelerate the shrinking speed, a suitable strategy for the
detection and analysis of clusters (connected components of ${\cal
A}_i$) is used together with the application of the restarting
criterion of [1] for choosing new centers for the disks in the
clusters by means of the computation of a convex hull.
 
\section{Outline of the algorithm}
The program tries to perform the most part of the computation in the
lowest working precision starting with the standard machine precision
and switching to multiple precision {\it only if needed}.  In this way,
for polynomials having only few ill conditioned roots the (more
expensive) mp version of the algorithm is applied just for the ill
conditioned roots with a substantial saving of the computation time.
This is made possible since our algorithm does use the ``clean'' 
information given by the coefficients of $p(x)$ at each step
of the computation, in fact, {\em it does not compute the explicit
factorization} of the polynomial.  Algorithms based on explicit
factorizations need a high precision computation of the coefficients
of the factors obtained just at the first splitting stage even though
the factors have only few ill conditioned roots.

This feature of using at each stage the information of the
coefficients of $p(x)$ is particularly useful for sparse polynomials
where the computation of $p(x)$ becomes very cheap.  In fact, we have
implemented a sparse version of the Ruffini-Horner rule for the
computation of $p(x)$ and of $p'(x)$. The higher speed of this sparse
version can be exploited at each level of our algorithm for polynomial
roots.

Moreover, the program allows also to use polynomials defined by
straight line programs. For instance the polynomial $p_k$ of degree
$n=2^k-1$ defined as
\begin{equation}\label{mandel}
\begin{array}{l}
p_0(x)=1 \\
p_{i}(x)=x(p_{i-1}(x))^2+1, \ \ i=1,2,\ldots, k,
\end{array}
\end{equation}
having roots strictly related to the Mandelbrot set, can be computed
with just $3k$ arithmetic operations and the computation is more
accurate.

Another feature of our algorithm is that it allows to use only those
digits of the input coefficients that are strictly needed for the
requested approximation. More specifically, it is typical in symbolic
computation to encounter polynomials with integer coefficients having
hundreds or even thousands digits. Very often, in order to separate or
to count the (real) roots of such polynomials only few leading digits
of the coefficients are enough.  MPSolve, starts the computation
with using few leading digits of the coefficients, and doubles the
number of them at each step of the computation until the requested
goal is reached.  This leads to a substantial saving of the time cost
in many concrete situations with respect to other algorithms based on
explicit factorization. Typically, for the latter class of algorithms,
the useful information that initially is
confined in the leading digits of the coefficients, is scrambled
 in all the (many) digits of the coefficients of the factors  after
just the first factorization step.

\section{Notes on the implementation}
The program is written in ANSI C, makes use of the GMP multiprecision
package, and of the extended types package MT (More Types) written by
G. Fiorentino.  Besides the type {\tt double}, and {\tt cplx}
(implementing complex numbers having standard C {\tt double} real and
imaginary parts) the package uses the {\tt rdpe} and the {\tt cdpe}
types consisting in real and complex numbers where the real components
are represented as a pair (mantissa, exponent), with a C {\tt double}
mantissa and a {\tt long int} exponent. The latter types are needed to
perform computation within the precision (and with a slightly higher
cost) of that of {\tt double}, but allow a much wider range that
avoids overflow and underflow situations. The notation ``dpe'', which
means Double Precision and Exponent, is borrowed by the MPFun package
by D. Bailey.  For higher precision computation MPSolve uses the
multiprecision types of GMP and the multiprecision complex type {\tt
mpc} implemented by G. Fiorentino.

In order to achieve a high speed of computation, almost all the
sub-algorithms of the package have been implemented in three slightly
different versions: the ``float'' version, the ``dpe'' version and the
``mp'' version.

\section{Using  MPSolve}
The executable program of the package MPSolve is called unisolve and can be
used as:
\begin{center}
 {\tt unisolve {\em options} {\em input\_file}}
\end{center}

If the {\tt \em input\_file} is missing then unisolve will read from 
the standard input stream (typically the keyboard).

\subsection*{Options and directives}
The program unisolve allows many {\tt \em options} to change
its default behavior,
here is the list:
\begin{description}
  \item[{\tt -G}{\em x}:] goal of the computation. Three goals are possible:
   \begin{description}
     \item{\tt i}: isolate (default)
     \item{\tt a}: approximate
     \item{\tt c}: count
   \end{description}
   For instance, {\tt unisolve -Gi } isolates the roots of the polynomial
   \item[{\tt -S}{\em x}:] set where the program looks for the roots.
       The available sets are
   \begin{description}
     \item{\tt a}: all the complex plane (default)
     \item{\tt l}: left half plane $\{z\in{\bf C}:\quad Re(z)<0\}$
     \item{\tt r}: right half plane $\{z\in{\bf C}:\quad Re(z)>0\}$
     \item{\tt u}: upper half plane $\{z\in{\bf C}:\quad Im(z)>0\}$
     \item{\tt d}: (down) lower half plane $\{z\in{\bf C}:\quad Im(z)<0\}$
     \item{\tt i}: inside the unit disk $\{z\in{\bf C}:\quad |z|<1\}$
     \item{\tt o}: outside the unit disk $\{z\in{\bf C}:\quad |z|>1\}$
     \item{\tt R}: Real line $\{z\in{\bf C}:\quad Im(z)=0\}$
     \item{\tt I}: Imaginary line $\{z\in{\bf C}:\quad Re(z)=0\}$
   \end{description}
    For instance,
   {\tt unisolve -Gc -SR } counts the real roots of the polynomial,
   {\tt unisolve -Gi -So } isolates the roots out of the unit disk.
   \item[{\tt -D}{\em x}:] the program detects if the root is real/imaginary.
    For a detected real root the imaginary part is written as `0',
    similarly, the real part of a detected imaginary root is written as `0'.
   The options are:
   \begin{description}
     \item{\tt r}: detect real
     \item{\tt i}: detect imaginary
     \item{\tt b}: detect both real and imaginary roots.
     \item{\tt n}: do not detect (default)
   \end{description}

   \item[{\tt -M}{\em x}:] the program detects if there are multiple  roots
   there are two options detect {\tt -M+}, do not detect {\tt -M- } (default).
   \item[{\tt -d}:] if this switch is set then the program writes, in the 
    standard error, information about the execution of the program.
    By using the modifier {\tt -d1} all output goes to the standard output
    stream.
   \item[{\tt -i}{\em x}:] input precision in digits.  
        If, say, the option {\tt -i100 } is set, then it is assumed
        that the coefficients of the input polynomial have 100 correct
        digits. This value overwrites the precision assigned in the
        input file. The option {\tt -i0 } (default) means infinite
        precision.

\item[{\tt -o}{\em x}:] output precision.
        If, say, the option {\tt -o100 } is set, then the program will
        isolate/approximate/count the roots of the polynomial up to
        the precision of 100 digits. The default value is 30 digits.
        The detection of multiplicities and of real/imaginary roots,
        that are based on the computation of the Mahler bound, are
        affected by this value.

\item[{\tt -O}{\em x}:] output format. The options for {\em x} are 
                (here re and im denote real and imaginary parts of
                the root, respectively):
        \begin{description}
           \item{\tt c}:  compact form (default): (re, im)       
           \item{\tt b}:  bare format: re $\backslash$tab  im 
           \item{\tt g}:  "Gnuplot", low precision format: re im, well
                           suited as gnuplot input
           \item{\tt v}:  verbose: Root(i) = re $\pm$ I im
           \item{\tt f}:  full information: (re, im), a posteriori 
                          error bound, status of the approximation.   
        \end{description}

With the option {\tt -Of} the real and imaginary parts are written
integrally with no filtering of the correct digits. The status of the
approximation is a triple of characters with the following meaning:
\begin{itemize}
  \item First character: (cluster, isolated, approximated, multiplicity)
   \begin{description}
     \item{m:}  multiple root    
     \item{i:}  isolated root 
     \item{a:}  approximated single root (relative error less than
                 $\epsilon_{out}$)
     \item{o:}  approximated cluster of roots (relative error less
                 than $\epsilon_{out}$) 
     \item{c:}  cluster of roots not yet approximated (relative error
                 greater than $\epsilon_{out}$) 
   \end{description}
  \item Second character: (real/imaginary roots)
     \begin{description}
       \item{R:} real root
       \item{r:} non-real root
       \item{I:} imaginary root
       \item{i:} non-imaginary root
       \item{w:} uncertain real/imaginary root
       \item{z:} non-real and non-imaginary root
     \end{description}
  \item Third character (set $\cal S$)
     \begin{description} 
       \item{i:} root in $\cal S$
       \item{o:} root out of  $\cal S$
       \item{u:} root  uncertain
      \end{description}
\end{itemize}


\item[{\tt -Lp}{\em x}:] maximum number of  iterations per packet (default 10).
     The Ehrlich-Aberth iterations are grouped into packets. The number of 
     iterations per packet can be modified.
 
\item[{\tt -Li}{\em x}:]  maximum number of global iterations (default 100). 
The number of packets of the Ehrlich-Aberth iteration can be modified.
Increasing this number may be needed for dealing with exceptionally hard 
polynomials.

\end{description}

\subsection*{The input file}
The {\tt \em input\_file} must  contain:
\begin{itemize}
\item some optional lines of comments, introduced by the character `{\tt !}'.
\item a string of three characters {\tt abc} where the first character
of the string specifies the way the polynomial is assigned, the second
and the third specify the type of the coefficients.
\item the input precision of the coefficients, in digits (0 means
infinite precision).
\item degree and coefficients of the polynomial coded according to the
  values of {\tt abc}.
\end{itemize}

\noindent
More precisely:
\begin{verbatim}
       a=`s' means sparse polynomial; 
       a=`d' means dense polynomial;
       a=`u' means polynomial provided by the user by 
             means of a c program in the file mps_usr.c 
             (a sample of this file is provided for the 
             Mandelbrot polynomial);

       b=`r' means real coefficients; 
       b=`c' means complex coefficients;
    
       c=`i' means integer real/imaginary parts;
       c=`q' means rational real/imaginary parts;
       c=`b' means bigfloat real/imaginary parts;
       c=`f' means float real/imaginary parts.
\end{verbatim}

\noindent
For {\tt a}=`s' the file contains:
\begin{verbatim}
       n: the degree of the polynomial;
       n_coeff: the number of the nonzero coefficients;
       i:    indices of the generic nonzero coefficient;
       a_i:  the coefficient of x^i;
       ...   the order of this list of coefficients is 
             not relevant.
\end{verbatim}

\noindent
For {\tt a}=`d' the file contains:
\begin{verbatim}
       n: the degree of the polynomial
       a_0 constant coefficient 
       a_1 coefficient of x
       ...
\end{verbatim}

\noindent
For {\tt b}=`r' the coefficients are stored as single number for
{\tt c} $\neq$ `q' while for {\tt c}=`q', they are stored as
\begin{verbatim}
       num: numerator
       den: denominator
\end{verbatim}

\noindent
For {\tt b}=`c' the coefficient are stored as pairs provided that
{\tt c} $\neq$ `q', i.e.,
\begin{verbatim}
       re: real_part
       im:  imaginary_part
\end{verbatim}

\noindent
otherwise, for {\tt c}=`q' they are stored as quadruples:
\begin{verbatim}
       real_numer: numerator of the real part 
       real_denom: denominator of the real part 
       imag_numer: numerator of the imaginary part
       imag_denom: denominator of the imaginary  part
\end{verbatim}

\noindent
For instance the following data defines the Kameny polynomial
\begin{verbatim}
! Kameny polynomial, sparse, complex, integer.
! p(x)=ic^3 x^7+ c^4 x^2- 6c^2x+9, c=10^6
sci
0
7
4
0
 9
 0
1
 -6000000000000
 0
2
 10000000000000000000000010
 0
7
 0
 1000000000000000000
\end{verbatim}

\noindent
The following data define the Wilkinson polynomial of degree 20.
\begin{verbatim}
! Wilkinson polynomial n=20: p(x)=\prod_{i=0}^n (x-i)
dri
0
20
2432902008176640000
-8752948036761600000
13803759753640704000
-12870931245150988800
8037811822645051776
-3599979517947607200
1206647803780373360
-311333643161390640
63030812099294896
-10142299865511450
1307535010540395
-135585182899530
11310276995381
-756111184500
40171771630
-1672280820
53327946
-1256850
20615
-210
1
\end{verbatim}

\noindent
Finally, the data
\begin{verbatim}
uri
0
511
\end{verbatim}

\noindent
define the user polynomial (by default it is the Mandelbrot polynomial)
of degree 511.

The program writes to the standard output stream, and, if the {\tt -d}
switch is set, also to the standard error stream.  The information written
in the standard output are: the approximations of the roots, where
each root is written as {\tt (real\_part, imaginary\_part)}.  The
information written to the standard error stream concern the
execution of the program, moreover, also the approximations of the
roots are reported, where the roots are not replaced by 0.0 if the
relative error is greater than 1.  Also the absolute error bound of
each approximation is reported. This information may be useful in the
case of polynomials with approximate input.

\section{The test suite}

Some input test polynomials are provided in the directory {\tt Data};
the directory {\tt Results} contains output files for some of the test
polynomials.

The test polynomials have the following features:
\begin{itemize}
\item The polynomials from {\tt spiral10} to {\tt spiral30} are defined by
\[\begin{array}{l}
p(x) = (x+1) (x+1+a) (x+1+a+a^2) ... (x+1+a+a^2+... +a^n), \\
 a = {\bf i}/1000, ~~   {\bf i}^2 = -1,~~   n = 10, 15, 20, 25, 30.
\end{array}
\] 
The polynomials from  {\tt geom1\_10} to {\tt geom1\_40} and from
 {\tt geom2\_10} to {\tt geom2\_40} are defined by
$p(x) = (x+1) (x+a) (x+a^2) ... (x+a^{n-1})$, 
with $a=100{\bf i}$ and $a={\bf i}/100$. respectively.
The polynomials from  {\tt geom3\_10} to {\tt geom3\_80} are defined by
$4^{n(n+1)/2}\prod_{i=1}^n(x-1/4^i)$, for $n=10,20,40,80$.
The polynomials from  {\tt geom4\_10} to {\tt geom4\_80} are defined by
$\prod_{i=1}^n(x-4^i)$ for $n=10,20,40,80$.

\item The polynomials from {\tt easy100} to {\tt easy6400} are defined by
$p(x)=\sum_{i=0}^n (i+1)x^i$, for $100, 200, 400, 800, 1600, 3200, 6400$. 

\item
The polynomials from {\tt mandel31} to {\tt mandel1023} are the Mandel
polynomials of degrees from 31 to 1023, defined in section 1 and
assigned in terms of their coefficients.
\item
The polynomials from {\tt umand31} to {\tt umand2047} are the Mandel
polynomials of degrees from 31 to 2047, assigned in terms of a
straight line program (file {\tt mps\_usr.c}).
\item
The polynomials from {\tt exp50} to  {\tt exp400} are the 
truncation of the
exponential series to the degree $n=50,100,200,400$, respectively, i.e.,
$\sum_{i=0}^n x^i/i!$.

\item Curz polynomials of degree 20 40 80 160 {\tt curz20}--{\tt curz160}
defined by:
$p_{j+1}(x)=x\sum_{i=0}^{j-1}p_{j-i}(-1)^i/(i+1) +(-1)^j/(j+1)$,
$j=1,2,\ldots$,
for $p_1=1$.

\item
The polynomials {\tt kam1\_1}, {\tt kam1\_2}, {\tt kam1\_3} belong to the
class
\[
(x-3c^2)^2+icx^7,
\]
for $c=10^{-6}, 10^{-20}, 10^{-70}$, respectively.  
The polynomials have two simple
complex roots with extremely small real separation (33, 113, 385,
common digits respectively) and imaginary parts close to zero (of
modulus about $10^{-44}, 10^{-149}, 10^{-524}$, respectively).  The
remaining roots are well separated.
\item
The polynomials {\tt kam2\_1}, {\tt kam2\_2}, {\tt kam2\_3} belong to the
class
\[
    (c^2x^2-3)^2+ic^2x^9,
\]
for $c=10^{6}, 10^{20}, 10^{70}$, respectively.  
The polynomials have four simple
complex roots with small imaginary parts ($10^{-27},
10^{-90},10^{-315}$, respectively,) in two close sets (21, 70, 247
common digits).
The remaining roots are well separated.

\item
The polynomials {\tt kam3\_1}, {\tt kam3\_2}, {\tt kam3\_2} belong to the
class
\[
    (c^2x^2-3)^2+c^2x^9,
\]
for $c=10^{6}, 10^{20}, 10^{70}$, respectively.
The polynomials have two extremely close real roots (20, 70, 247
common digits) plus a complex pair with extremely small imaginary
parts \\
($10^{-27}, 10^{-90}, 10^{-315}$)
The remaining roots are well separated, one of them is real.

\item
The polynomial {\tt kam4} is defined as follows.
\[
   x^{14}+2 \cdot 10^{24}\cdot x^{11}+10^{48}\cdot x^8+4 \cdot x^7-4
   \cdot10^{24}\cdot x^4+4
\]
The polynomial has a cluster of 8 roots of modulus $\tt 1.89\cdot
10^{-6}$, a cluster of 6 roots having modulus $\tt 10^8$. The cluster
made up by 8 roots has: two real roots having 21 common digits
two pairs of complex roots having small real parts, and imaginary
parts with 21 common digits
The cluster made up by 6 roots has two very close real roots matching
in their first 28 digits, and 2 very close complex pairs whose real
and imaginary parts match in their first 27 digits.

\item The polynomials from {\tt mig1\_20}, to {\tt mig1\_500},
and from {\tt mig1\_50\_1} to\break  {\tt mig1\_500\_1} are Mignotte-like
polynomials in the class                  
\[
    x^n+(a x+1)^m, \ \ n>>m>1, \ |a|>1,
\]
These polynomials have $m$ very close roots around $-1/a$. 
The polynomials have been chosen by setting $m=3$ and $m=31$. 


\item
The polynomial {\tt mig2} is defined by
\[
(x^{n-k}+(a x+1)^m)(a x+1)^k, \ \ n>>m>1, \ n>>k>1,\ |a|>1,
\]
for $ n=20, m=3, k=3, a=100 i$. 
The polynomial has a cluster of $m$ very close roots around $-1/a$
where $-1/a$ is also a multiple root of multiplicity $k$.

\item The polynomials {\tt lar1} and {\tt lar1\_200} are defined as 
\[
x^{n}+10^{300}x^{14}+x^5+1, \ \ 
\]
have a coefficient with modulus close to the largest number
representable as double in the IEEE standard.  Overflow
conditions/failure are likely for programs that do not use extended
arithmetic.  The polynomials are obtained with $a=1$ and $n=20,200$,
respectively.

\item The polynomial {\tt lar2} is defined by
\[
  x^{n}+x^{11}+10^{300}x+10^{-300},  n=20.
\]
The polynomial has coefficients that have modulus close to the largest
and to the smallest numbers representable as double in the IEEE
standard.  Overflow conditions/failure are likely for programs that do
not use extended arithmetic.  The polynomial, even if its coefficients
are representable as double, has a root which is too small in order to
be represented as a nonzero double.
\item the polynomial {\tt lar3} is defined by
\[
 10^{-200} x^{n}+10^{200} x^{n-1}+10^{200},  \ \  n=20
\]

The polynomial cannot be normalized in the standard IEEE arithmetic
without generating an overflow condition.  The polynomial, even if its
coefficients are representable as double, has a root which cannot be
represented as a double without generating an overflow condition.
Failure is likely for programs that do not use extended arithmetic.

\item The polynomial {\tt lar4} is defined by 
$
10^{20}+10^{2000} x^{19} +
10^{-1600} x^{23}
$.

\item The polynomial {\tt lar5} is defined by 
$
10^{2000} 
+
10^{2000} x^{19}
+
10 x^{20}
$.

\item Polynomial {\tt lsr1}
\[
\begin{array}{l}
1 + 30000x + 300000000x^2 + 1000000000000x^3 + x^{200} + 400000000
  x^{298} + \\ 12000000040000 x^{299} + 120000001200000001 x^{300} +
  \\ 400000012000000030000 x^{301} + 40000000300000000 x^{302} + \\
  1000000000000 x^{303} + 400000000 x^{498} + \\ 40000 x^{499 }+
  x^{500}
\end{array}
\]
The polynomial has a cluster of two very close roots of large moduli
and a cluster of three very close roots of small moduli.

\item Polynomial {\tt lsr2} given by
\[
 x^{500} +100^{200} x^{499}+ 10^{500} x^{495}+ 10^{499} x^{480}+10^3
 x^3+3\cdot 10^2 x^2+30 x+1
\]
The polynomial has roots with large and small moduli.

\item Polynomial {\tt lsr3} defined by
$1+x*30+x^2*300+x^3*1000+x^{480}*10^{495}
+x^{495}* 10^{500}
+x^{499}*10^{200}
+x^{500}$.


\item Polynomials {\tt lsr4\_1},{\tt lsr4\_2}, {\tt lsr4\_3} are given by  
$(x^{50}+1)(x^2+ax+a^{-1})$, $a=10^{10}, 10^{20}, 10^{40}$, respectively.


\item Polynomial {\tt lsr\_200}
$(x^{12}-(10^{20} x-1)^4)(1+(10^{20}+x)^4 x^8)(x^{200}+1)$. \\
 There are 4 roots of modulus $10^{-20}$ matching in the
  first 60 digits and 4 roots of modulus $10^{20}$ matching in the
  first 60 digits.


\item Polynomial {\tt lsr\_24}
 $(x^{12}-(10^{20} x-1)^4)(1+(10^{20}+x)^4 x^8)$.\\
  There are 4 roots of modulus $10^{-20}$ matching in the
 first 60 digits and 4 roots of modulus $10^{20}$ matching in the
 first 60 digits.

\item The polynomials from {\tt nroots50} to {\tt nroots1600} are of
the kind $z^n-1$, $n=50,100,200,400,800,1600$.

\item The polynomials from  {\tt wilk20} to {\tt wilk320} are the
  Wilkinson polynomial of degree 20,40,80,160,320, defined by
\[
   \prod_{i=1}^n(x-i), \ \ \ (n=20).
\]
Coefficients are very large and the roots are very ill-conditioned
already for $n=20$.  The condition number of the roots, for $n=20$
ranges from 420 to $10^{14}$.

\item The polynomial {\tt wilk\_mod1} is a modified Wilkinson 
polynomial defined as follows
\[
   ((x^{10}/10^{20}+(x-10)^2) \prod_{i=1}^{20}(x-i).
\]
The polynomial is obtained by multiplying Wilkinson's polynomial of
degree 20 with a Mignotte-like polynomial having two roots close to
10.  Besides the ill conditioning of the integer roots, there is a
cluster of two roots close to 10. The condition number of the root 10
is about $10^{26}$.

\item The polynomial {\tt wilk\_mod2} is a modified Wilkinson polynomial
defined as
\[
   \prod_{i=1}^{20}(x-i)(x-20)^2.
\]
Besides the ill conditioning of the roots, there is a multiple root.


\item The polynomials {\tt kir1\_10},  {\tt kir1\_20},  {\tt kir1\_40},
proposed by Peter Kirrinnis, are defined as follow
$ (z^4-1/16)^n(z^4-(1/2+eps)^4)16^n/eps^4$, $eps=1/4096$,
$n=10,20,40$ and have degree $4n+4$. They have multiple roots with
high multiplicity and some of them clustered.


\item The polynomials {\tt kir1\_10\_mod},  {\tt kir1\_20\_mod},  
{\tt kir1\_40\_mod},
are obtained by adding $z^{4n+4}$ to the Kirrinnis polynomials. In this
way they
are defined as follow
$ (z^4-1/16)^n(z^4-(1/2+eps)^4)16^n/eps^4+z^{4n+4}$, $eps=1/4096$,
$n=10,20,40$ and have degree $4n+4$. They have simple clustered roots.

\item The polynomial {\tt nektarios}, provided by Nektarios Psycaris,
has degree 648 and arises in the study of certain problems in physics. 

\item The polynomial {\tt trv\_m}, provided by Carlo Traverso, arises
from
the symbolic processing of a system of polynomial equations. It has
multiple roots.

\item The polynomial {\tt katsura8} arises from the symbolic
preprocessing
of a system of polynomial equations. 

\item The polynomials {\tt chrm*} are chromatic polynomials.
\end{itemize}


The HTML file {\tt bench.htm} in the {\tt Doc} directory reports the timings,
in seconds, of the programs NSolve, fsolve, polroots and unisolve of
the packages Mathematica, Maple, Pari and MPSolve respectively. All
computations have been performed on a Pentium 200 MMX with 64 Mb
of physical RAM.

We performed computations with several goals, like isolating all the
roots, approximating all of them with 50 or 1000 digits, counting the
real roots or those in the unit disk or those in the right half plane,
approximating only real roots, etc.

Failures, denote with an `F' in the table, occurred in Maple for
internal errors like `{\em unable to compute the norm}', or for
unsuccessful computation of the roots. In Pari we had failures only
due to stack overflows, while in Mathematica we had unsuccessful
computation of the roots.  We point out that the Mathematica command
NSolve, when used with the $d$ digits of working precision, say {\tt
NSolve[p[x]==0,x,100]} for $d=100$, does not guarantee $d$ correct
output digits. In fact the number of correct output digits depends on
the conditioning of the root. In particular, in the case of a root of
multiplicity $m$, in order to have at least $k$ correct digits we had
to choose $d=km$, i.e., to type {\tt NSolve[p[x]==0, x, k*m]}.

Timings are complete only for the package MPSolve.

\section{User defined polynomials.}
The file {\tt mps\_user.c} contains the ``user defined polynomial'',
more precisely, the programs {\tt fnewton\_usr, dnewton\_usr,
mnewton\_usr} that are the float, dpe, mp versions, respectively of
the same computation.  Actually, the float version is not used by {\tt
unisolve} (indeed, the float version cannot be as robust as the dpe
version for overflow reasons) and has been inserted as a more readable
example.

Given on input the value of $x$, each program in the file {\tt
mps\_user.c} must compute:
\begin{itemize}
\item The value of $p(x)$
\item The value of $p'(x)$
\item An upper bound $\sigma$ to the number $|p(x)-fl(p(x))|$, 
      where $fl(\cdot)$ denotes the floating point computation of
      the expression between parentheses,
\end{itemize}
and output the following values:
\begin{itemize}
\item the Newton correction {\tt corr}$=p(x)/p'(x)$;
\item the inclusion radius {\tt rad}$=n(|p(x)|+\sigma)/|p'(x)|$;
\item the Boolean value {\tt again}$= (|fl(p(x))| > \sigma)$, where
      if {\tt again} is true then the Newton iteration is
      continued. Observe that {\tt again} true means that the
      computed value is affected by a relative roundoff error less
      than 1, i.e.,
      $|fl(p(x))-p(x)|/|fl(p(x))|<\sigma/fl(p(x))<1$. Therefore it
      brings still information.
\end{itemize} 

Indeed, computing $p(x)$ and $p'(x)$ is almost straightforward if the
formal relations defining the polynomial are known.  For instance, for
the Mandelbrot polynomial defined in (\ref{mandel}), it is immediate
to write a program computing $p(x)$. For the computation of $p'(x)$ we
may use the following relation obtained by differentiating
(\ref{mandel}).

\[
\begin{array}{l}
p'_0=0 \\ p'_{i}=p_{i-1}^2+2\cdot x\cdot p_{i-1}\cdot p'_{i-1}
\end{array}
\]
In order to implement the stop condition and to give a guaranteed
a-posteriori error bound, a rounding error analysis of the computation
must be performed.

Let us show this for the Mandelbrot polynomial.  Let $\widetilde
p=fl(p)$ denote the value obtained by applying (\ref{mandel}) in
floating point arithmetic with machine precision $\mu$ and recall the
following relations (compare [1]) $fl(a*b)=a*b(1+\epsilon)$ and
$fl(a+b)=(a+b)(1+\delta)$, where $|\epsilon|<2\sqrt 2\mu$ and
$|\delta|<\mu$, which hold for complex values $a,b$ that do not
generate overflow/underflow in the computation of $a+b$ and $a*b$
(here $\mu$ is the machine precision).

By applying the above relations to (\ref{mandel}) we easily obtain
\[
\begin{array}{l}
\widetilde p_{i}=x\widetilde p_{i-1}^2(1+\theta_i)+(1+\sigma_i) \\
|\theta_i|<6.7\mu +O(\mu^2),\ |\sigma_i|<\mu+O(\mu^2)
\end{array}
\]
Subtracting both the members of the above relation and of
(\ref{mandel}), ignoring $O(\mu^2)$ terms yields the following
expressions that allow us to compute a bound to the absolute error
$e_i=|p_i-\widetilde p_i|$
\[
\begin{array}{l}
e_0=0\\ e_{i}<2|x|\cdot e_{i-1}|\widetilde
p_{i-1}|+7.7(|x|\cdot|\widetilde p_{i-1}|+1)\mu
\end{array}
\]
The rounding error analysis can be performed in a rough way obtaining
very weak bounds, or in a very accurate way obtaining sharp bounds.
 
The more accurate is the roundoff error bound, the more efficient is
the stop criterion.

\section*{References}

[1] D. A. Bini and G. Fiorentino, Design, Analysis, and Implementation
of a Multiprecision Polynomial Rootfinder, Numerical Algorithms,
23, 2000, 127--173.

[2] D. A. Bini, Numerical computation of polynomial zeros by means of
Aberth's method, Numerical Algorithms, 13, 1996, 179--200.

\noindent
[3] D. A. Bini and G. Fiorentino, A multiprecision implementation of a
poly-algorithm for univariate polynomial zeros, Proc. of The POSSO
Workshop on Software, Paris, 1995, J.C. Faug\`ere, J. Marchand, R.
Rioboo editors.

\noindent
[4] F. Rouillier, RUR, FRISCO Report 3.3.2.3.1.

\end{document}