Solving univariate polynomials




1  Aberth method

D. Bini & G. Fiorentino
This is a connection to the C-implementation by D. Bini and G. Fiorentino for solving univariate polynomials. See [Bini97abrt] for more details.

upoly/solve/Aberth.H

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)*x30+(-465)*x29+(103385)*x28+(-14631225)*x27+(1480321269)*x26+(-114009431445)*x25
+(6949189247325)*x24+(-344092707928125)*x23+(14097793282984515)*x22+(-484338676679532675)*x21+(14090257524223082475)*x20+(-349600545868057540875)*x19
+(7435941626111727234855)*x18+(-136055808711963322871175)*x17+(2145883249334501452139775)*x16+(-29197210605623737977801375)*x15+(342563613932937660652700640)*x14+(-3460266110493898677911394000)*x13
+(30006513636556697864066736800)*x12+(-222457423246962063058403076000)*x11+(1401937624086807501691142239744)*x10+(-7454161471690660700139655157760)*x9+(33114629767614997850763390570240)*x8+(-121365366674745136523074652102400)*x7
+(360930788158836812805614538878976)*x6+(-851899888505423112503184251412480)*x5+(1547794975254719737111781253120000)*x4+(-2070792202024594683660866641920000)*x3+(1902893785240928209998216560640000)*x2+(-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.25877017821792895931348326448096581026e-61)
#2: (0.19999999999999999999999999999999999984e1, -0.46556680990471469891214099942031669987e-76)
#3: (0.40000000000000000000000000000000364749e1, -0.20535888801067914054504324258151119759e-62)
#4: (0.29999999999999999999999999999999998682e1, -0.1690014420565147123865015402332071313e-67)
#5: (0.60000000000000000000000000010308179885e1, -0.17069344456114038582494063531601900797e-55)
#6: (0.49999999999999999999999999999761733257e1, 0.56117336544939643771944799678521116096e-47)
#7: (0.99999999999999999999999784826963231212e1, 0.1677037890900643849669321014094637549e-45)
#8: (0.90000000000000000000000032519985952991e1, -0.10509095875125750363385727455212235666e-47)
#9: (0.14000000000000000000008402553702753828e2, 0.2015435089101158787266077152350015994e-45)
#10: (0.69999999999999999999999999991601077172e1, 0.25390980354637881208074621189882365307e-37)
#11: (0.10999999999999999999999885325082034558e2, -0.57235100952987539063160947223702612529e-40)
#12: (0.79999999999999999999999999431433217699e1, -0.19467375500487617082017239633722323927e-41)
#13: (0.12999999999999999999996347819253591034e2, 0.36353917880566144751326965843923149935e-48)
#14: (0.12000000000000000000000970795413911516e2, 0.10465433766150597832348156110497136095e-28)
#15: (0.16000000000000000000121692969865195616e2, 0.61282400956055226193665133857607039139e-44)
#16: (0.16999999999999999999732193650094651459e2, 0.2309095887500468055812263958690829855e-38)
#17: (0.17999999999999999999682383096594987337e2, 0.16812616094710630162783072442274842416e-30)
#18: (0.15000000000000000000002189951417373643e2, -0.7490939520500240247661306054365703158e-32)
#19: (0.19999999999999999999296586876141711257e2, 0.85176035248872305472653559861172840824e-31)
#20: (0.18999999999999999999529788281178441923e2, 0.53282144207137213154959172264956273312e-40)
#21: (0.21999999999999999999297348006978163749e2, 0.36943424133849245652733336879058978001e-37)
#22: (0.22999999999999999999287914975831207769e2, 0.29732741191732769710887783935183503301e-41)
#23: (0.25999999999999999999832045321006631775e2, 0.20580343010841541263617393462696565414e-34)
#24: (0.20999999999999999998631601379682898985e2, 0.85416495423790731130368144063784408372e-38)
#25: (0.28999999999999999999999511858249512079e2, 0.58775478356019121465377999997238523646e-29)
#26: (0.25000000000000000000182522016583431122e2, 0.1702818085845692419455710925989157286e-36)
#27: (0.2800000000000000000001097792243810599e2, -0.27922288574837944256539061438507685528e-31)
#28: (0.24000000000000000000512743725846821381e2, -0.18359428755560974665678223574819808583e-30)
#29: (0.27000000000000000000036044721387544108e2, 0.49609681736270557227520758047099445496e-31)
#30: (0.29999999999999999999999913709517042125e2, 0.82974904262406279750260379614458531495e-31)

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 fn, we compute the successive squares f2n.
An iterative process based on the Sebastiao e Silva method

  1. Choose a starting value f of A.
  2. Iterate fn+1 = (fn)2/ ||(fn)2)|| in A.
If there exists a simple root z maximizing |f|, we obtain
limn-> infty fn+1=ez.
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 ft = (1,t,t2,..., td-1), that is the coefficient vector in the Horner basis of p(x) - p(t)/ (x-t)p'(t). We have ft(z) = -p(t)/ z-t for all the roots of p, so that the evaluation of ft at the roots of p is maximal when |z-t| is minimal. If we want the smallest root of p, we start the process from f0. 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)/ x-z. The second deflation process, called implicit deflation, is as follows. One root z is computed with a starting value f0, and the method produced a first idempotent ez. To get a second root, we choose any g of A (possibly f again) and start the method from the initial value (g-gez). Since (g-g  ez)(z) = 0, z is no longer dominant, and the process cannot converge again to ez. Thus we will obtain another idempotent ez' 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 h-h(ez+ez'). 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

upoly/solve/Sebastio.H

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 B1,p(x,x0) 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 x0, using Sebastiao e Silva method. If there are two roots closest to x0, 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 x0. 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:
n
S
i=0
xi
i!
See the examples supplied by D.A. Bini in the FRISCO test suite [PolTS].
  1. The polynomial is obtained by truncating the series of ex.
  2. 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
  1. Polynomials with real random coefficients, of degrees 10000 and 100000.
  2. Polynomials with complex random coefficients, of degrees 10000, 100000 and 1000000.
  3. Exponential polynomials Sixi/ i!, for some degrees of 100, 150, 200.
  d x0 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, x0 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 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
Df(x0,t)= |f(x0)| -
d
S
k=1
|f(k)(x0)|
k!
tk
where f is a polynomial of degree d. We can check that if Df(x0,r)³ 0, the interval [x0-r,x0+r] does not contain a root of f.

upoly/solve/Exclusion.H

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.