- Introduction
- Finding bounds on the roots
- First Cauchy theorem
- Second Cauchy theorem
- Third Cauchy theorem
- Lagrange-MacLaurin theorem
- Laguerre method
- Laguerre second method
- Newton method
- Newton theorem
- Joyal bounds
- Pellet method
- Global implementation
- Kantorovitch theorem

- Bounds on the product and sum of roots

- Maximum number of real roots
- Number of real roots

- Separation between the roots

- Analyzing the real roots
- Analyzing the real part of the roots
- Utilities

Analyzing univariate polynomials

Finding bounds on the roots

First Cauchy theorem

Mathematical background

with . Let

Then the modulus of any root of verify:

int Cauchy_First_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`Bound`: the interval on the absolute value of the roots

int Cauchy_First_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound);This procedure returns an interval such that for any polynomial in the set and for any root of this polynomial . A failure code of 0 is returned if

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; Num=Cauchy_First_Bound_Interval(3,Coeff,Bound); Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=Cauchy_First_Bound_Interval(3,Coeff_App,Bound);In the first case the procedure find that the absolute value of the roots lie in [0.6,4] while in the second case the range is [0.58,4.44444].

Second Cauchy theorem

Mathematical background

with . Let with the strictly negative coefficients of . Then all the positive real roots of verify [13]:

Note that if all the roots are negative according to Descartes Lemma (see section 5.5.1).

int Cauchy_Second_Bound_Interval(int Degree,VECTOR &Coeff,double *Bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`bound`: the upper bound of the positive roots

int Cauchy_Second_Bound_Inverse_Interval(int Degree,VECTOR &Coeff,double *Bound);This procedure fail and returns 0 if

int Cauchy_Second_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL Bound);This procedure will return a failure code of 0 if

int Cauchy_Second_Bound_Inverse_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL Bound);In that case if

A global procedure enable to get at the same time the upper and lower bound of the positive roots:

int Cauchy_Second_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound); int Cauchy_Second_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Lower,INTERVAL &Upper);In the latter case the interval lower bound is in

int Cauchy_Second_Bound_Negative_Interval(int Degree,VECTOR &Coeff,double *Bound); int Cauchy_Second_Bound_Negative_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound);while the upper bound of the negative roots can be determined using:

int Cauchy_Second_Bound_Negative_Inverse_Interval(int Degree,VECTOR &Coeff,double *Bound); int Cauchy_Second_Bound_Negative_Inverse_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Bound);Both the upper and lower bound of the negative roots can be found using:

int Cauchy_Second_Bound_Negative_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound); int Cauchy_Second_Bound_Negative_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Lower,INTERVAL &Upper);

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; Num=Cauchy_Second_Bound_Interval(3,Coeff,&bound); Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=Cauchy_Second_Bound_Interval(3,Coeff_App,Bound);In the first case we find that all the roots are lower than 2 and in the second that the roots have bounds [2.44444,2.44444]. If we have used the lower bound procedures we will have found that all the roots are positive.

Third Cauchy theorem

Mathematical background

int Cauchy_Third_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound)

int Cauchy_All_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Bound)will return the result of a successive application of the first and third Cauchy bounds.

Lagrange-MacLaurin theorem

Assume and let () be the first negative coefficients of (if has no negative coefficients then there is no positive real root).

The upper bound of the value of the positive real root is:

where is the greatest absolute value of the negative coefficients of ,[3],[13].

If we define:

Then the upper bound of the positive real roots of is the lower bound of the positive real root of . Consequently if and are computed for the polynomial then

int MacLaurin_Bound_Interval(int Degree,VECTOR &Coeff,double *Bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`bound`: the upper bound of the positive roots

On success the return code is 1. There is also a procedure for the interval polynomial:

int MacLaurin_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL Bound);This procedure fail and returns 0 if

It is also possible to determine the lower bound of the positive roots using the procedures:

int MacLaurin_Bound_Inverse_Interval(int Degree,VECTOR &Coeff,double *Bound); int MacLaurin_Bound_Inverse_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound);In the later case if

int MacLaurin_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound) int MacLaurin_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Lower, INTERVAL &Upper)In the latter case the interval lower bound is in

It is also possible to determine the lower bound of the negative roots using the procedures:

int MacLaurin_Bound_Negative_Interval(int Degree,VECTOR &Coeff,double *Bound); int MacLaurin_Bound_Negative_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound);Similarly it is possible to determine the upper bound of the negative roots using the procedures:

int MacLaurin_Bound_Negative_Inverse_Interval(int Degree,VECTOR &Coeff,double *Bound); int MacLaurin_Bound_Negative_Inverse_Interval(int Degree, INTERVAL_VECTOR &Coeff,INTERVAL &Bound);To get directly both lower and upper bound of the negative roots you may use:

int MacLaurin_Bound_Negative_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound) int MacLaurin_Bound_Negative_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Lower, INTERVAL &Upper)

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; Num=MacLaurin_Bound_Interval(3,Coeff,&bound); Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=MacLaurin_Bound_Interval(3,Coeff_App,Bound);In the first case we find that all the roots are lower than 4 and in the second that the roots have bounds [4.44444,4.44444]. If we have used the lower bound procedures we will have found that all the roots are greater than -1.

Laguerre method

with , . Let's define the sequence :

If it exists a real such that for all in [0,], then for all . Consequently all the roots of are lower than [13]. To find the following scheme can be used:

- let be such that
- let the smallest integer such that either or and
- if then substitute by such that and go to 2
- return

int Laguerre_Bound_Interval(int Degree,VECTOR &Coeff,double sens,int MaxIter,double *bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`sens`: a positive real indicating the increase of in the scheme`MaxIter`: the maximum number of iteration. If this number is exceeded the procedure returns 0.`bound`: the upper value of the real root

The lower bound of the root may be determined by:

int Laguerre_Bound_Inverse_Interval(int Degree,VECTOR &Coeff1,double amp_sens, int MaxIter,double *bound);We have also a procedure which determine upper and lower bound for the real roots:

int Laguerre_Bound_Interval(int Degree,VECTOR &Coeff,double sens,int MaxIter,INTERVAL &Bound);All the real roots lie within

int Laguerre_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,double sens, int MaxIter,INTERVAL &Bound);This procedure fail and returns 0 if

int Laguerre_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,double sens, int MaxIter,INTERVAL &Lower,INTERVAL &Upper);In that case

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; Num=Laguerre_Bound_Interval(3,Coeff,0.1,50,Bound); Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=Laguerre_Bound_Interval(3,Coeff_App,0.1,50,Lower,Upper);leads to

Laguerre second method

Assume that has real roots. Then [23] all roots are contained in the interval whose end-point are given by the two solutions of the quadratic equation:

int Laguerre_Second_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`Bound`: upper and lower bound on the roots

int Laguerre_Second_Bound_Interval(int Degree, INTERVAL_VECTOR &Coeff,INTERVAL &Lower,INTERVAL &Upper);In that case

Newton method

with , . Newton theorem state that if it exists such that and for all the derivative of , with , then is an upper bound of the positive roots of . To find the following scheme can be used:

- let be such that
- let the smallest integer such that either or and
- if then substitute by such that and go to 2
- return

In the implementation we define as the smallest real which satisfy
. If we find such that
, then we
increase by a given positive value `sens` and start again. We
limit the number of iteration of the scheme by giving a maximal value
for the number of iteration:

int Newton_Bound_Interval(int Degree,VECTOR &Coeff1,double amp_sens,int MaxIter,double *bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`sens`: a positive real indicating the increase of in the scheme`MaxIter`: the maximum number of iteration. If this number is exceeded the procedure returns 0.`bound`: the upper value of the real root

int Newton_Bound_Inverse_Interval(int Degree,VECTOR &Coeff1,double amp_sens, int MaxIter,double *bound);We have also a procedure which determine upper and lower bound for the real roots:

int Newton_Bound_Interval(int Degree,VECTOR &Coeff1,double amp_sens, int MaxIter,INTERVAL &Bound);All the real roots lie within

int Newton_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,double sens, int MaxIter,INTERVAL &Bound);This procedure fail and returns 0 if

int Newton_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,double sens, int MaxIter,INTERVAL &Lower,INTERVAL &Upper);In that case

Newton theorem

with real roots. The roots of are such that [4]:

int Newton_Second_Bound_Interval(int Degree,VECTOR &Coeff,double *bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`bound`: upper bound on the absolute value of the roots

int Newton_Second_Bound_Inverse_Interval(int Degree,VECTOR &Coeff,double *bound);An upper and lower bound on the absolute value of the roots may be compute by:

int Newton_Second_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound);In that case if

Similar procedures exist for interval polynomial:

int Newton_Second_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound); int Newton_Second_Bound_Inverse_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Bound); int Newton_Second_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &L, INTERVAL &U);where

`Bound`: upper or lower bound on the absolute value of the root of the polynomials`L`,`U`: if`L`=[a,b] then the absolute value of the roots of any polynomial in the set is greater than a while some polynomial have root greater than b. Conversely if`U`=[a,b] the absolute value of the roots of any polynomial in the set is lower than b while some polynomial have root lower than a.

Joyal bounds

int Joyal_Bound_Interval(int Degre,INTERVAL_VECTOR &Coeff,INTERVAL &Bound)returns in

Pellet method

Now assume that we have calculated the coefficients of the polynomial where are known quantities. If we determine a polynomial that has 2 positive roots in the range , then has roots in the disk . Hence the absolute value of the real part of the roots is bounded by . As we get that the real part of the root satisfies . As this shows that has roots whose real part is greater than .

int Pellet(int Degre,INTERVAL_VECTOR &Coeff,double R)where

int Global_Positive_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound); int Global_Negative_Bound_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`Bound`: the bound on the positive or negative roots

- if a bound have been previously determined, then the step
size is fixed to except if the global variable
`Step_Laguerre_Interval`has been defined to be a double not equal to 0, in which case this variable is used as the step size. - if a bound as not been determined the step size is fixed to 1
except if the global variable
`Step_Laguerre_Interval`has been defined to be a double not equal to 0, in which case this variable is used as the step size. - if the step if lower than the global variable
`Min_Step_Laguerre_Interval`(which is 0.1 by default) the step is substituted by`Min_Step_Laguerre_Interval`

int Global_Positive_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Lower,INTERVAL &Upper); int Global_Negative_Bound_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL &Lower,INTERVAL &Upper);where

Both procedures for real roots have been regrouped in

void ALIAS_Find_Bound_Polynom(int Degree, INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), INTERVAL_VECTOR &PALL,INTERVAL &Space, INTERVAL &Bound)where

`PALL`: the first component is the bound for the real roots the polynomial while the following are the value of the parameters for which the bounds will be obtained`Space`: a default bound for the real roots; may be used to indicate that we are looking only for bounds on the positive or negative real roots by specifying a positive or a negative lower bound`Bound`: bounds for the real roots

Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=Global_Positive_Bound_Interval(3,Coeff_App,Lower,Upper);which leads to

The above procedure may give some sharp bounds. For example consider the Wilkinson polynomial of order 22 (which has as roots 1,2,,21,22) we get 0 negative roots while the positive roots are bounded by [0.790447,22.1087].

The test program `Test_Bound_UP` enable to test the bound
procedures for any polynomial. This program take as first argument
the name of a file giving the
coefficients of the polynomial by increasing power of the
variable. The program will print the bounds determined by all the
previous procedures and then the bounds determined using the global
implementation. Then the same treatment will be applied on the
interval polynomial whose coefficients are intervals centered at the
coefficients find in the file with a width of 0.2

Kantorovitch theorem

int Kantorovitch(int Degree,VECTOR &Coeff,REAL Input,double *eps)with

`Degree`: degree of the polynomial`Coeff`: coefficients of the polynomial ordered along increasing degree`Input`: center of the interval in which an unique solution occur`eps`: half width of the solution interval

int Kantorovitch(int Degree,INTERVAL_VECTOR &Coeff,REAL Input,double *eps)If this procedure returns 1, then any polynomial in the set of interval polynomial has an unique solution in the interval [

int Kantorovitch_Fast_Safe(int Degree,INTERVAL_VECTOR &Coeff,REAL Input,double *eps)in which "safe" interval value of the coefficients have been pre-computed.

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; Num=Kantorovitch(3,Coeff,1.2,&eps); Coeff_App(1)= INTERVAL(-3.1,-2.9);Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9);Coeff_App(4)=INTERVAL(0.9,1.1); Num=Kantorovitch(3,Coeff_App,1.2,&eps);In the first case

Newton relations

and the real and complex roots of . Let define as

We have [4]:

This procedure enable to compute the elements . The syntax is:

VECTOR SumN_Polynomial_Interval(int Degree,VECTOR &Coeff)with:

`Degree`: degree of the polynomial`Coeff`: coefficients of the polynomial ordered along increasing degree

INTERVAL_VECTOR SumN_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff)which returns intervals including the . This procedure returns 0 if 0 is included in the leading interval.

Viète relations

and the real and complex roots of . Let define as:

We have [4]:

This procedure enable to compute the elements ( is the sum of the roots, the product of the roots). The syntax is:

VECTOR ProdN_Polynomial_Interval(int Degree,VECTOR &Coeff)with:

`Degree`: degree of the polynomial`Coeff`: coefficients of the polynomial ordered along increasing degree

INTERVAL_VECTOR ProdN_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff)which returns intervals including the . This procedure returns 0 if 0 is included in the leading interval.

Descartes Lemma

with . Let the sequence and the the number of change of sign in this sequence. Then the number of positive real roots, counted with their order of multiplicity is equal to with [13].

INT Descartes_Lemma_Interval(int Degree,VECTOR &Coeff)with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree

There is an implementation of this method for interval polynomial. Here it is necessary to introduce an additional parameter to indicate the confidence we have in the result. The procedure is implemented as:

INT Descartes_Lemma_Interval(int Degree,INTERVAL_VECTOR &Coeff,int *Confidence);

- 1: the result is exact, for all polynomials in the set the number of positive real roots is
- : the number of positive real roots for all the polynomials in the set is

INT Descartes_Lemma_Negative_Interval(int Degree,VECTOR &Coeff) INT Descartes_Lemma_Negative_Interval(int Degree,INTERVAL_VECTOR &Coeff,int *Confidence)

Budan-Fourier method

and its n-th derivative. Let the interval the interval in which we are looking for roots. We assume that , and . We construct the sequence from which we exclude the 0 element. Similarly we construct the sequence (a special treatment has to be applied for the zero element of , see [13]). Let the number of change of sign in and the number of change of sign in . Then the number of real roots of in , counted with their order of multiplicity, is or lower than this number by an even number.

INT Budan_Fourier_Interval(int Degree,VECTOR &Coeff,INTERVAL In) INT Budan_Fourier_Interval(int Degree,INTEGER_VECTOR &Coeff,INTERVAL In)with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree, which may be`REAL`or`INT``In`: the interval in which we are looking for the number of roots

- -1:
- -2: , the polynomial may be factored
- -3:

INT Budan_Fourier_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL In,int *Confidence)where

- 1: the result is exact, so if is the return code of the algorithm the number of reals roots of the interval polynomial is with
- : if is the return code of the algorithm the number of reals roots of the interval polynomial is

INT Budan_Fourier_Safe_Interval(int Degree,VECTOR &Coeff, INTERVAL In,INTERVAL &NbRoot);The procedure returns 1 in case of success and an interval for the number of roots. If

INT Budan_Fourier_Fast_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff, INTERVAL In,INTERVAL &NbRoot);Another safe procedure is:

INT Budan_Fourier_Interval(int Degree,INTEGER_VECTOR &Coeff,int Inf,int Sup)where the coefficients and the bounds are integers.

Coeff(1)= -3;Coeff(2)=2;Coeff(3)=-1;Coeff(4)=1; P=INTERVAL(0.0,2.0); Num= Budan_Fourier_Interval(3,Coeff,P);

Coeff_App(1)= INTERVAL(-3.1,-2.9); Coeff_App(2)=INTERVAL(1.9,2.1); Coeff_App(3)=INTERVAL(-1.1,-0.9); Coeff_App(4)=INTERVAL(0.9,1.1); Num= Budan_Fourier_Interval(Degree,Coeff_App,P,&Confidence);

Sturm method

Let be the number of change of sign in the first sequence and be the number of change of sign in the second sequence. Then the number of distinct real roots of in the interval is [14] if . Note that a multiple roots count only for one root with this method.

The drawback of Sturm method is that the absolute value of the coefficients increase quickly when computing the sequence. Numerical rounding errors may then affect the result. The alternate method of Budan-Fourier (see section 5.5.2) is less sensitive to rounding errors although it provides less information than Sturm method.

This procedure determines the number of real roots of an univariate polynomial in a given interval. It is implemented as:

int Sturm_Interval(int Degree,VECTOR &Coeff,INTERVAL &In); int Sturm_Interval(int Degree,INTEGER_VECTOR &Coeff,INTERVAL &In);

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree, either real or integers`In`: the interval in which we are looking for the number of roots

int Sturm_Safe_Interval(int Degree,VECTOR &Coeff,INTERVAL &In,INTERVAL &NbRoot);which returns 1 in case of success and an interval

int Sturm_Interval(int Degree,INTEGER_VECTOR &Coeff,int Inf,int Sup);This procedure will return -1 if all the coefficients are 0 and -2 if at some point of the process an integer larger than the largest machine integer is encountered (

and we are looking for roots in the interval [0.5,0.7] the procedure

The test program `Test_Nb_Root_Up` enable to test Budan-Fourier
and Sturm method on any polynomial whose coefficients are defined in a
file, by increasing degree of the power of the unknown.

Du Gua-Huat-Euler theorem

If all the roots of are real then [13]:

Conversely if this relation does not hold for some (for example if a coefficient is 0 while its neighbor have same sign), then there are complex roots.

int Huat_Polynomial_Interval(int Degree,VECTOR &Coeff);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree

int Huat_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff);This procedure will return 1 if all the polynomials in the set have their roots real, -1 if some polynomials in the set have all roots real but some others have complex roots and 0 if all the polynomials in the set have complex roots.

Separation between the roots

Rump theorem

and let be:

The minimal distance between the real roots of is such that [26]:

int Min_Sep_Root_Interval(int Degree,VECTOR &Coeff,double &min);with

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients of the polynomial in increasing degree`min`: the lower bound on the distance between two real roots

int Max_Sep_Root_Interval(int Degree,VECTOR &Coeff,double &max);while upper and lower bounds may be determined with:

int Bound_Sep_Root_Interval(int Degree,VECTOR &Coeff,INTERVAL &Bound);There is also a procedure to determine a lower bound for interval polynomial:

int Min_Sep_Root_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL &Lower);If

int ALIAS_Min_Max_Is_Root(int Degree, int NbParameter, int HasInterval, INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), int Iteration, INTERVAL_VECTOR &Par, double *Root, int Type, int (* Solve_Poly)(double *, int *,double *))may be used to rest if a polynomial in a set may have as real root one of two pre-defined value.

`Degree`: the degree of the polynomial`NbParameter`: the number of parameters that appear in the coefficient of the polynomial`HasInterval`: 1 if the coefficient include intervals not defined by parameters, 0 otherwise`TheCoeff`: a procedure to calculate the coefficients of the polynomial being given range for the parameters`Iteration`: the maximum number of box that may be used by the algorithm`Par`: the ranges for the parameters`Root`: we look for polynomial whose real part of the root is either Root[0] or Root[1]`Type`:- -1 : if in a box we have found a polynomial with a root lower than Root[0]; we look for a polynomial in the box whose real root is exactly Root[0]. If no such polynomial is found the box is eliminated
- 1 : if in a box we have found a polynomial with a root greater than Root[1]; we look for a polynomial in the box whose real root is exactly Root[1]. If no such polynomial is found the box
- 2: if in a box we have found a polynomial with a root greater than Root[1] and a polynomial with a real root greater than Root[0]; we look for a polynomial in the box whose real root is exactly Root[1] and a polynomial whose real root is exactly Root[0]. If no such polynomials are found the box is eliminated
- 0: general case;

`Solve_Poly`: a procedure to solve polynomial with double floating point coefficients. The first argument are the coefficients, the second the degree of the polynomial and the third the real roots

For some applications it may be interesting to determine if a polynomial in a set defined by a parametric polynomial has the real part of one of its roots equal to a pre-defined value. This may be done by using the procedure

int ALIAS_Is_Root_RealPart(int Degree, int NbParameter, int HasInterval, INTERVAL_VECTOR (* TheCoeff)(INTERVAL_VECTOR &), INTERVAL_VECTOR (* TheCoeffCentered)(INTERVAL_VECTOR&,double a), int Iteration, INTERVAL_VECTOR &Par, double *Root, INTERVAL_VECTOR (* EvaluateComplex)(int,int,INTERVAL_VECTOR &), int (* Simp)(INTERVAL_VECTOR &))where

`Degree`: the degree of the polynomial`NbParameter`: the number of parameters that appear in the coefficient of the polynomial`HasInterval`: 1 if the coefficient include intervals not defined by parameters, 0 otherwise`TheCoeff`: a procedure to calculate the coefficients of the polynomial being given range for the parameters`TheCoeffCentered`: a procedure to calculate the coefficients of the polynomial`Iteration`: the maximum number of box that may be used by the algorithm`Par`: the ranges for the parameters`Root`: we look for polynomial whose real part of the root is either Root[0] or Root[1]`EvaluateComplex`: a procedure that returns 4 intervals. The input interval vector has dimension`NbParameter+1`elements, the first one being the parameters, the last one being . The procedure should evaluate and should return the real part of , the complex part of , the real part of and the complex part of`Simp`: an optional simplification procedure that returns -1 if the polynomials in a set cannot roots with real part equal to Root[0] or Root[1]

Addition of two polynomials

The following procedures enable to add polynomials with real,interval or integer coefficients:

VECTOR Add_Polynomial_Interval(int n1,VECTOR &Coeff1,int n2,VECTOR &Coeff2); INTERVAL_VECTOR Add_Polynomial_Interval(int n1,INTERVAL_VECTOR &Coeff1, int n2,INTERVAL_VECTOR &Coeff2); INTEGER_VECTOR Add_Polynomial_Interval(int n1,INTEGER_VECTOR &Coeff1, int n2,INTEGER_VECTOR &Coeff2);with:

`ni`: degree of the polynomial`Coeffi`: the coefficients of the polynomials in increasing degree

INTERVAL_VECTOR Add_Polynomial_Safe_Interval(int n1,VECTOR &Coeff1,int n2,VECTOR &Coeff2);

Multiplication of two polynomials

The following procedures returns the degree of the product of two polynomial either with real or integer coefficients:

int Degree_Product_Polynomial_Interval(int n1,VECTOR &Coeff1,int n2,VECTOR &Coeff2); int Degree_Product_Polynomial_Interval(int n1,INTEGER_VECTOR &Coeff1,int n2, INTEGER_VECTOR &Coeff2);Then you may use the following procedure to compute the product of polynomial either with real,integer or interval coefficients:

VECTOR Multiply_Polynomial_Interval(int n1,VECTOR &Coeff1,int n2,VECTOR &Coeff2); INTEGER_VECTOR Multiply_Polynomial_Interval(int n1,INTEGER_VECTOR &Coeff1, int n2,INTEGER_VECTOR &Coeff2); INTERVAL_VECTOR Multiply_Polynomial_Interval(int n1,INTERVAL_VECTOR &Coeff1, int n2,INTERVAL_VECTOR &Coeff2);All these procedures return the coefficients of the product, the leading term being 0 only if one polynomial is equal to 0.

To take into account the rounding errors you may use:

INTERVAL_VECTOR Multiply_Polynomial_Safe_Interval(int n1,VECTOR &Coeff1,int n2,VECTOR &Coeff2);

Evaluation of a polynomial

The evaluation of a polynomial for a given value is implemented as:

REAL Evaluate_Polynomial_Interval(int Degree,VECTOR &Coeff,REAL P) INTERVAL Evaluate_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff,REAL P) INTERVAL Evaluate_Polynomial_Interval(int Degree,VECTOR &Coeff,INTERVAL P) INTERVAL Evaluate_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL P) REAL Evaluate_Polynomial_Interval(int Degree,INTEGER_VECTOR &Coeff,REAL P); int Evaluate_Polynomial_Interval(int Degree,INTEGER_VECTOR &Coeff,INT P);with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients (which can be`REAL`,`INT`or`INTERVAL`of the polynomial in increasing degree`P`: the point at which we want to compute the polynomial. It may be`REAL`,`INT`or`INTERVAL`.

- the value of a polynomial with
`REAL`coefficients at a point - the interval value of a polynomial with interval coefficients at a real point
- the interval value of a polynomial with
`REAL`coefficients at for an interval value of the unknown. The first and second order derivative of the polynomial are used to get sharp bounds - the interval value of a polynomial with interval coefficients for an interval value of the unknown. The first and second order derivative of the polynomial are used to get sharp bounds

Evaluation in centered form

[11977379510.705810546875,11977410100.791748046875]A better evaluation may be obtained if we use a centered form of the polynomial. Consider the polynomials:

where is a real. For an appropriate choice of the we may have for all . The procedures:

VECTOR Coeff_Polynomial_Centered_Interval(int Degree,VECTOR &Coeff,REAL P); INTERVAL_VECTOR Coeff_Polynomial_Centered_Interval(int Degree, INTERVAL_VECTOR &Coeff,REAL P);enable to compute the for the centered form at

INTERVAL_VECTOR Coeff_Polynomial_Centered_Safe_Interval(int Degree,VECTOR &Coeff,REAL P); INTERVAL_VECTOR Coeff_Polynomial_Centered_Fast_Safe_Interval(int Degree, INTERVAL_VECTOR &Coeff,REAL P);which return safe value for the coefficients (in the second form we assume that you have pre-computed safe value for the coefficients of the polynomial using the procedure described in section 5.9.10).

Then we may use the procedures:

REAL Evaluate_Polynomial_Centered_Interval(int Degree,VECTOR &Coeff,REAL Center,REAL P); INTERVAL Evaluate_Polynomial_Centered_Interval(int Degree,VECTOR &Coeff,INTERVAL P); INTERVAL Evaluate_Polynomial_Centered_Interval(int Degree,INTERVAL_VECTOR &Coeff, REAL Center,REAL P); INTERVAL Evaluate_Polynomial_Centered_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL P);These procedures return the evaluation of the polynomial at

Safe evaluation of a polynomial

INTERVAL Evaluate_Polynomial_Safe_Interval(int Degree,VECTOR &Coeff,REAL P); INTERVAL Evaluate_Polynomial_Safe_Interval(int Degree,VECTOR &Coeff,INTERVAL P); INTERVAL Evaluate_Polynomial_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff,REAL P); INTERVAL Evaluate_Polynomial_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL P);

For example for the Wilkinson polynomial at order 15 the
evaluation for 15.1 leads to the interval:

If "safe" intervals have been pre-computed (for example by using the procedure described in section 5.9.10) for the coefficients you may use:

INTERVAL Evaluate_Polynomial_Fast_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff,REAL P); INTERVAL Evaluate_Polynomial_Fast_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL P);which are faster. If you want to use a centered form you may use:

INTERVAL Evaluate_Polynomial_Centered_Safe_Interval(int Degree,VECTOR &Coeff,REAL P); INTERVAL Evaluate_Polynomial_Centered_Safe_Interval(int Degree,VECTOR &Coeff,INTERVAL P); INTERVAL Evaluate_Polynomial_Centered_Fast_Safe_Interval(int Degree, INTERVAL_VECTOR &Coeff,REAL P); INTERVAL Evaluate_Polynomial_Centered_Fast_Safe_Interval(int Degree, INTERVAL_VECTOR &Coeff,INTERVAL P);in which the two last forms assume that you have pre-computed a safe value of the coefficients of the polynomial.

The evaluation of a polynomial with real coefficients at a complex point may be performed with:

INTERVAL_VECTOR Evaluate_Complex_Poly(int deg, INTERVAL_VECTOR &Coeff, INTERVAL &XR,INTERVAL &XI)where

Sign of a polynomial

int Sign_Polynomial_Interval(int Degree,VECTOR &Coeff,REAL P) int Sign_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff,REAL P) int Sign_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff,REAL P) int Sign_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff,INTERVAL P)with:

`Degree`: degree of the polynomial`Coeff`: the`Degree+1`coefficients (which can be`REAL`or`INTERVAL`of the polynomial in increasing degree`P`: the point or interval at which we want to compute the polynomial

This procedure returns 1 if the polynomial is positive, -1 if it is
negative, 0 if it is 0 and 2 if the sign is either affected by
numerical errors (if `Coeff` and `P` are
`REAL`) or is not constant (if
`Coeff` or `P` are `INTERVAL_VECTOR`).

Derivative of a polynomial

VECTOR Derivative_Polynomial_Interval(int Degree,VECTOR &Coeff) INTEGER_VECTOR Derivative_Polynomial_Interval(int Degree,INTEGER_VECTOR &Coeff); VECTOR Nth_Derivative_Polynomial_Interval(int Degree,VECTOR &Coeff,int n)These procedures enable to get the coefficients of the first and

Due to rounding errors there may be errors in the coefficients provided by the previous procedures. The procedure:

INTERVAL_VECTOR Derivative_Polynomial_Safe_Interval(int Degree,VECTOR &Coeff);will return an

INTERVAL_VECTOR Derivative_Polynomial_Fast_Safe_Interval(int Degree,INTERVAL_VECTOR &Coeff);

If we deal with polynomial with interval coefficients the following procedures may be used:

INTERVAL_VECTOR Derivative_Polynomial_Interval(int Degree,INTERVAL_VECTOR &Coeff) INTERVAL_VECTOR Nth_Derivative_Polynomial_Interval(int Degree, INTERVAL_VECTOR &Coeff,int n)

Euclidian division

The procedure:

void ALIAS_Euclidian_Division(int Degree,INTERVAL_VECTOR &Coeff1, INTERVAL_VECTOR &Coeff2,INTERVAL &Residu,INTERVAL &a)computes the division of the polynomial

Another procedure that perform the same process if `a` is a double
is

void Divide_Single(int Degree,INTERVAL_VECTOR &Coeff1,double a, INTERVAL_VECTOR &Coeff2, INTERVAL &Residu)For dividing a polynomial by the product of polynomials of type so that

and evaluating the polynomial under this form for some interval for you may use

INTERVAL Divide_Evaluate_Multiple(int Degree,INTERVAL_VECTOR &Coeff1, int n,double *a,INTERVAL_VECTOR &X)where

`n`: the number of terms`a`: the value of

For the Euclidian division by an arbitrary polynomial you may use

int Divide_Polynom(int Degree,INTERVAL_VECTOR &Coeff,int DegreeDivisor, INTERVAL_VECTOR &CoeffDivisor, INTERVAL_VECTOR &CoeffQuo, INTERVAL_VECTOR &CoeffRem)The polynomial of degree

Note that `Divide_Polynom` and `Divide_Evaluate_Multiple`
may be used to perform the deflation of an univariate polynomial as
soon as approximate roots for the polynomial have been found.
This may be useful to decrease the computation time for solving a
polynomial or for numerically instable polynomial such as the
Wilkinson polynomial. A combination of Rouche filtering and of
deflation allows to solve Wilkinson polynomial of order 18, while the
general solving procedure will fail starting at order 13.

A special Euclidian division is the division of a polynomial by its derivative. This can be done with

void Quotient_UP_Derivative(int Degree,VECTOR &Coeff,VECTOR &CoeffD,VECTOR &Quo, VECTOR &Rem)where has coefficients

Expansion of

The coefficients of the expansion of the polynomial may be obtained with the procedure

INTERVAL_VECTOR Power_Polynomial(INTERVAL &a,int n)which returns the coefficients

Let which may also be written as where is some fixed or interval value. The coefficient may be calculated with the procedure

int Derive_Polynomial_Expansion(int n,INTERVAL_VECTOR &Coeff,INTERVAL &a,INTERVAL_VECTOR &B)which returns 1 if the have been successfully calculated (0 otherwise). The coefficients of are

Unitary polynomial

Let be a polynomial and be `maxroot` the maximal modulus of
the root of . From we may derive a polynomial such that the
roots of have a modulus lower or equal to 1 and if is a root
of then `maxroot` is a root of .
The calculation of the coefficients of may be done with the
procedure

int Unit_Polynom(int Degree,INTERVAL_VECTOR &Coeff, double maxroot,INTERVAL_VECTOR &CoeffU)where

Safe evaluation of a vector

INTERVAL_VECTOR Evaluate_Coeff_Safe_Interval(int n,VECTOR &Coeff);with:

`n`: size of the vector`Coeff`: the vector

2018-07-25