Solving univariate polynomials 
1 Aberth method
D. Bini & G. Fiorentino
This is a connection to the Cimplementation by D. Bini and
G. Fiorentino for solving univariate polynomials.
See [Bini97abrt] for more details.
template< class T>
VAL< VectStd< array1d< Scl< MPC> > > *> Solve(const T &t, const Aberth & A)
Solve the univariate polynomial of type T, by the Aberth method.
The complex roots are output in a vector of Scl< MPC> complex numbers.
The default constructor of the class Aberth has to be used when we assume
no error on the input coefficients and we want the default precision (2*DBL_DIG digits) on the output roots.
The constructor Aberth(i) is used when we assume an input precision of i
digits and the default output precision for the approximation of the
roots.
The constructor Aberth(i,o) is used when we assume an input precision of i
digits and we want an output precision of o digits on the roots.
Example:
#include <string.h>
#include "upoly.H"
#include "bignum.H"
#include "util/Clock.H"
#include "upoly/solve.H"
typedef Scl<MPZ> C;
typedef Scl<MPC> CC;
typedef UPolyDense<upar<C> > upoly;
typedef VectStd<array1d<CC> > Vect;
int main (){
upoly q; cin >>q; cout<<q<<endl;
(1)*x^{30}+(465)*x^{29}+(103385)*x^{28}+(14631225)*x^{27}+(1480321269)*x^{26}+(114009431445)*x^{25} 
+(6949189247325)*x^{24}+(344092707928125)*x^{23}+(14097793282984515)*x^{22}+(484338676679532675)*x^{21}+(14090257524223082475)*x^{20}+(349600545868057540875)*x^{19} 
+(7435941626111727234855)*x^{18}+(136055808711963322871175)*x^{17}+(2145883249334501452139775)*x^{16}+(29197210605623737977801375)*x^{15}+(342563613932937660652700640)*x^{14}+(3460266110493898677911394000)*x^{13} 
+(30006513636556697864066736800)*x^{12}+(222457423246962063058403076000)*x^{11}+(1401937624086807501691142239744)*x^{10}+(7454161471690660700139655157760)*x^{9}+(33114629767614997850763390570240)*x^{8}+(121365366674745136523074652102400)*x^{7} 
+(360930788158836812805614538878976)*x^{6}+(851899888505423112503184251412480)*x^{5}+(1547794975254719737111781253120000)*x^{4}+(2070792202024594683660866641920000)*x^{3}+(1902893785240928209998216560640000)*x^{2}+(1059681761389533859949327155200000)*x 
+(265252859812191058636308480000000)

Clock c; c.start();
Vect r =Solve(q, Aberth());
cout<<"Roots:"<<endl;
for (int i = 0; i < q.degree(); i++) cout<<"#"<<i+1<<": "<<r[i]<<endl;
}
Roots:
#1: (0.1e1, 0.25877017821792895931348326448096581026e61)
#2: (0.19999999999999999999999999999999999984e1, 0.46556680990471469891214099942031669987e76)
#3: (0.40000000000000000000000000000000364749e1, 0.20535888801067914054504324258151119759e62)
#4: (0.29999999999999999999999999999999998682e1, 0.1690014420565147123865015402332071313e67)
#5: (0.60000000000000000000000000010308179885e1, 0.17069344456114038582494063531601900797e55)
#6: (0.49999999999999999999999999999761733257e1, 0.56117336544939643771944799678521116096e47)
#7: (0.99999999999999999999999784826963231212e1, 0.1677037890900643849669321014094637549e45)
#8: (0.90000000000000000000000032519985952991e1, 0.10509095875125750363385727455212235666e47)
#9: (0.14000000000000000000008402553702753828e2, 0.2015435089101158787266077152350015994e45)
#10: (0.69999999999999999999999999991601077172e1, 0.25390980354637881208074621189882365307e37)
#11: (0.10999999999999999999999885325082034558e2, 0.57235100952987539063160947223702612529e40)
#12: (0.79999999999999999999999999431433217699e1, 0.19467375500487617082017239633722323927e41)
#13: (0.12999999999999999999996347819253591034e2, 0.36353917880566144751326965843923149935e48)
#14: (0.12000000000000000000000970795413911516e2, 0.10465433766150597832348156110497136095e28)
#15: (0.16000000000000000000121692969865195616e2, 0.61282400956055226193665133857607039139e44)
#16: (0.16999999999999999999732193650094651459e2, 0.2309095887500468055812263958690829855e38)
#17: (0.17999999999999999999682383096594987337e2, 0.16812616094710630162783072442274842416e30)
#18: (0.15000000000000000000002189951417373643e2, 0.7490939520500240247661306054365703158e32)
#19: (0.19999999999999999999296586876141711257e2, 0.85176035248872305472653559861172840824e31)
#20: (0.18999999999999999999529788281178441923e2, 0.53282144207137213154959172264956273312e40)
#21: (0.21999999999999999999297348006978163749e2, 0.36943424133849245652733336879058978001e37)
#22: (0.22999999999999999999287914975831207769e2, 0.29732741191732769710887783935183503301e41)
#23: (0.25999999999999999999832045321006631775e2, 0.20580343010841541263617393462696565414e34)
#24: (0.20999999999999999998631601379682898985e2, 0.85416495423790731130368144063784408372e38)
#25: (0.28999999999999999999999511858249512079e2, 0.58775478356019121465377999997238523646e29)
#26: (0.25000000000000000000182522016583431122e2, 0.1702818085845692419455710925989157286e36)
#27: (0.2800000000000000000001097792243810599e2, 0.27922288574837944256539061438507685528e31)
#28: (0.24000000000000000000512743725846821381e2, 0.18359428755560974665678223574819808583e30)
#29: (0.27000000000000000000036044721387544108e2, 0.49609681736270557227520758047099445496e31)
#30: (0.29999999999999999999999913709517042125e2, 0.82974904262406279750260379614458531495e31)
2 Sebastio a Sylva method
The Sebastiao e Silva method ([Card95]) is an iterative process
for computing a selected approximate root of a polynomial p. It is
essentially the same as the power method of matrices, but it is
accelerated: instead of computing the successive powers f^{n}, we
compute the successive squares f^{2n}.
An iterative process based on the Sebastiao e Silva
method
 Choose a starting value f of A.
 Iterate f_{n+1} = (f_{n})^{2}/ (f_{n})^{2}) in A.
If there exists a simple root z maximizing f, we obtain
lim_{n> infty} f_{n+1}=e_{z}.
Here, we have choosen the euclidian norm for the vectors. The process stops
if the error between two iteration is less than a value e,
fixed in advance.
Moreover, if the process is written in the Horner basis, an additionnal
interest is to have a simple way to select the root which is the nearest to a
given point, as illustred below.
Let us assume that, for a given point t
in K, there is only one simple root z closest to t. We start
the algorithm from the vector f_{t} = (1,t,t^{2},..., t^{d1}), that
is the coefficient vector in the Horner basis of p(x) 
p(t)/ (xt)p^{'}(t).
We have f_{t}(z) = p(t)/ zt for all the roots of p,
so that the evaluation of f_{t} at the roots of p is maximal when
zt is minimal. If we want the smallest root of p, we start
the process from f_{0}. If we want the largest we should reverse p and
look for its smallest root (see [Card95]).
If we want to obtain other roots after computing a first root z, we
apply a ``deflation'' of the polynomial p. There are two different
deflation processes, the first one is the most natural; it is called the
explicit deflation, and consists in computing explicitly
p(x)/ xz. The second deflation process, called implicit
deflation, is as follows. One root
z is computed with a starting value f_{0}, and the method produced a
first idempotent e_{z}. To get a second root, we choose any g of
A (possibly f again) and start the method from the initial value
(gge_{z}). Since (gg e_{z})(z) = 0, z is no
longer dominant, and the process cannot converge again to e_{z}. Thus
we will obtain another idempotent e_{z'} where
z' is the dominant root associated to g, z will be ignored. For
the next step, we choose h (possibly f or g) and start the method from
the initial value hh(e_{z}+e_{z'}). Continuing in this way, we
will ultimately collect as many roots of p as desired (for more details see
[Card95]).
For this work we have choosen the implicit deflation.
2.1 Implementation
struct Sebastiao
Class which gives the parameters of the Sebastiao method:
nbroot number of roots that should be computed,
nbiter maximal number of iterations allowed,
epsilon tolerance to stop the iteration.
{
int nbroot; //number of roots.
unsigned nbiter; //number of iterations.
double epsilon; //tolerance.
Sebastiao():nbroot(1),nbiter(40),epsilon(Tolerance){}
Sebastiao(const Sebastiao & m):
nbroot(m.nbroot),nbiter(m.nbiter),epsilon(m.epsilon){}
Sebastiao(int n):
nbroot(n),nbiter(40),epsilon(Tolerance) {}
Sebastiao(int n, unsigned nit):
nbroot(n),nbiter(nit),epsilon(Tolerance) {}
Sebastiao(int n, unsigned nit, double eps):
nbroot(n),nbiter(nit),epsilon(eps) {}
};
template < class R> inline
VAL< UPolyQuot< R> *> StartingValue(UPolyDense< R> p,
const typename R::value_type x0)
Returns B_{1,p}(x,x_{0}) in the dual basis. Used to get a starting value
of the iterative process.
template< class R>
VAL< UPolyQuot< R> *> IterFrom(const UPolyQuot< R> & p0, Sebastiao mth)
Apply Sebastiao e Silva iteration from the starting polynomial p0.
It stops either if the norm of the distance between two iterations
is smaller than mth.epsilon, or if the number of iterations goes
up to mth.nbiter.
template < class R> inline
typename R::value_type SolveOne(UPolyDense< R> & p,
const typename R::value_type & x0, Sebastiao mth)
Compute the root of the polynomial p which is the closest (in modulus)
to x_{0}, using Sebastiao e Silva method. If there are two roots closest
to x_{0}, the sequence will have two accumulation points and the output
solution will not be correct. This is not detected for the moment.
template < class R> inline
typename R::value_type* Solve(UPolyDense< R> & p,
const typename R::value_type & x0, Sebastiao mth)
Compute k=mth.nbroot roots of the polynomial p, using Sebastiao
e Silva method. It starts with the root the closest to x_{0}.
If k=1, all the roots of p are computed.
2.2 Tests:
Here are some results of experimentation for the Sebastiao e Silva
method. Some polynomials with random coefficients real and complex have been
generated, and we applied the method for the computation of one root for a
degree of 10000 and 100000. We have also computed all the root for a
polynomial with
random complex coefficients of degree 10000. One can observe that, for a
polynomial p with random coefficients, as its degree increases, the roots
of p will come closer and closer to the unit circle of the complex plane.
One root has also been computed for some univariate polynomials called truncated exponential, and this for different degrees. These polynomials
are of the form:
See the examples supplied by D.A. Bini in the
FRISCO test suite [PolTS].
 The polynomial is obtained by truncating the series of e^{x}.

The roots of this polynomial are located along a curve. To have a chance to
obtain a correct result for a degree > 20, it is necessary to use a multiple
precision library. Here it is the class Scl<MPF> based on
GNU MP ([GMP]).
In the following table, we present the results for
 Polynomials with real random coefficients, of degrees 10000 and 100000.
 Polynomials with complex random coefficients, of degrees 10000, 100000
and 1000000.
 Exponential polynomials S_{i}x^{i}/ i!, for
some degrees of 100, 150, 200.

d 
x_{0} 
k 
T 
double 
10000 
0.0 
8 
22.03s 
double 
100000 
0.0 
8 
346.52s 
complex 
10000 
(0.4,0.6) 
6 
16.86s 
complex 
100000 
(0.4,0.6) 
6 
292.02s 
complex 
1000000 
(0.4,0.6) 
7 
2520s 
exp 
100 
(0.0,0.8) 
6 
45s 
exp 
150 
(0.0,0.8) 
6 
59s 
exp 
200 
(0.0,0.8) 
6 
101s 
where d is the degree of the polynomial, x_{0} the starting value for the
iterative process, k is the number of iterations to compute one root and
T is the total time of the computation. These experimentations have been
performed on a Dec workstation with 496M of local memory.
3 Exclusion methods
The methods of this family are based on a test which answers
 No root in the interval.
 I don't know.
If we are in the first case, we remove the interval. Otherwise, the
interval is splitted into two subintervals, to which we apply the test
recursively.
The test implemented below use the function
D_{f}(x_{0},t)= f(x_{0})  


t^{k}

where f is a polynomial of degree d. We can check that
if D_{f}(x_{0},r)³ 0, the interval
[x_{0}r,x_{0}+r] does not contain a root of f.
template< class R>
typename R::value_type Exclude(const R & p,
const typename R::value_type & x0,
const typename R::value_type & r)
Test if the box B(x0,r) does not contain a (real) root.
If the result is positive, this is the case. Otherwise, we cannot say.
p is a univariate polynomial, x0 is the center of the box, and r
its radius
template< class POL, class C>
VAL< vector< C> *> Locate(const POL& p,const C& a,const C& b,const C & epsilon)
Output a sequence of points, such that the real roots
of the univariate polynomial p in the interval [a,b], are at a distance
less than epsilon, from one of these points.
template< class POL, class C> inline
VAL< vector< C> *> Locate(const POL& p, const box< C> & b, const C & epsilon)
The same but with a box b.