Stef Graillat
email: stef.graillat@lip6.fr

Fabienne Jézéquel
email: fabienne.jezequel@lip6.fr

Queiros Enzo

Maxime Martins

Spyropoulos

Computing multiple roots of polynomials in stochastic arithmetic with Newton method and approximate GCD

Keywords: approximate GCD, Discrete Stochastic Arithmetic, floating-point arithmetic, Newton method, numerical validation, rounding errors

Introduction

Polynomials appear in almost all areas of scientific computing. Generally the problem involved is to find roots of univariate polynomials. The wide range of use of polynomials needs to have fast and reliable methods to solve them. Roughly speaking, there are two general approaches: symbolic and numeric. The symbolic approach is based either on the theory of Gröbner basis or on the theory of resultants. For the numeric approach, we use iterative methods like Newton method or homotopy continuation methods. The symbolic algorithms give an exact representation of the result but can sometimes be slow. The numerical algorithms (performed in finite precision) are in general faster but give an approximate result. Our aim is to use numerical methods to get fast results without sacrificing accuracy.

If formulas are well-known for polynomials of degree less than 5, iterative methods are needed for greater degrees. In particular, Newton method is very used for its quadratic convergence. Other methods exist like Ehrlich-Aberth method, Durand-Kerner method or Laguerre method [START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF], [START_REF] Mcnamee | Numerical methods for roots of polynomials[END_REF]. In this article, we will only focus on Newton method.

When working with floating-point arithmetic, one major problem is to deal with ill-posed problems. Computing roots with multiplicities is an example of such a problem. Indeed, a small perturbation in the coefficients of a polynomial can change a double root into two single roots like for P ε (x) = (x -1) 2 -ε. If ε = 0 then P ε has 1 as a double root whereas for ε > 0, P ε has two single roots 1 ± √ ε. A similar phenomenon appears in the computation of the GCD of two polynomials. If P ε (x) = (x -1) 2 -ε and Q(x) = (x -1), then gcd(P ε (x), Q(x)) = x -1 for ε = 0 but gcd(P ε (x), Q(x)) = 1 for ε = 0. This has led to the development of different notions of approximate GCD or quasi-GCD [START_REF] Schönhage | Quasi-gcd computations[END_REF]- [START_REF] Nagasaka | Toward the best algorithm for approximate GCD of univariate polynomials[END_REF]. There are various definitions for approximate GCD but they closely follow the same concept. Given two polynomials P (x) and Q(x) the approximate GCD is somehow related to the computation of the exact GCD of polynomials P (x) and Q(x) close to P (x) and Q(x) with the objective that this GCD has the greatest possible degree.

Computing such approximate GCD is generally expensive in terms of execution time. In order to avoid this, we propose in this article to compute a stochastic GCD (denoted as stgcd) that makes it possible to smooth the computation of an approximate GCD using the classic Euclidean algorithm. It can be viewed as a kind of regularization technique that replaces floating-point arithmetic by stochastic arithmetic [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF].

This stochastic GCD will be used for deflation in Newton method. Given a polynomial P (x) with possible multiple roots, we will compute G(x) = st-gcd(P (x), P (x)). We will then use G(x) to compute the roots of P (x)/G(x) which normally has only single roots.

The contributions of the paper can be sum up as follows.

•

We use stochastic arithmetic to define a stochastic GCD that takes into account rounding errors when computations are performed in finite precision.

•

We use the new GCD to deflate polynomials in order to obtain polynomials with single roots.

•

We apply this to Newton's iterations and we show that we are able to accurately compute multiple roots of polynomials in finite precision.

The rest of the paper is organized as follows. In Section 2, we present how Discrete Stochastic Arithmetic (DSA) can be used to estimate rounding errors. Section 3 is devoted to the computation of multiple roots with DSA and the stochastic GCD. Numerical experiments are described in Section 4. Finally, concluding remarks are given in Section 5.

Estimation of rounding errors using Discrete Stochastic Arithmetic

This section presents the CESTAC method and Discrete Stochastic Arithmetic (DSA) as well as their implementations: CADNA and SAM. Those are well-known concepts and we heavily rely on [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF]. Further information can be found in [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF] and [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF].

Principles and validity of the CESTAC method

The CESTAC method [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF] enables one to estimate the rounding error propagation which occurs with floating-point arithmetic. This probabilistic method uses a random rounding mode: at each elementary operation, the result is rounded up (towards +∞) or down (towards -∞) with the probability of 50 %. Hence, with this random rounding mode, the same program run several times provides different results. Therefore, the computer's deterministic arithmetic is replaced by a stochastic arithmetic, where each arithmetic operation is performed N times before the next one is executed. The CESTAC method supplies us with N samples R 1 , . . . , R N of the computed result R. The value of the computed result R is then chosen as the mean value of {R i } and, if no overflow occurs, its number of correct digits (i.e. its number of digits not affected by rounding errors) can be estimated as

C R = log 10 √ N R στ β ( 1
)
where

R = 1 N N i=1 R i and σ 2 = 1 N -1 N i=1 R i -R 2 .
τ β is the value of Student's distribution for N -1 degrees of freedom and a confidence level 1 -β. In practice, β = 5% and N = 3. Therefore the number of correct digits is estimated within a 1 -β = 95% confidence interval. It has been shown [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF] that N = 3 is in some reasonable sense the optimal sample size. The estimation with N = 3 is more reliable than with N = 2 and increasing N does not significantly improve the quality of the estimation.

If both operands in a multiplication or the divisor in a division have no correct digits, the validity of C R is compromised [START_REF] Vignes | A stochastic arithmetic for reliable scientific computation[END_REF]. Therefore, the CESTAC method requires, during the execution of the user code, a dynamical control of multiplications and divisions, which is a so-called self-validation of the method.

Principles of DSA

The self-validation of the CESTAC method requires its synchronous implementation. Indeed, to enable the estimation of the accuracy, the samples which represent a result must be computed simultaneously. Discrete Stochastic Arithmetic (DSA) [START_REF] Vignes | Discrete Stochastic Arithmetic for validating results of numerical software[END_REF] has been defined from the synchronous implementation of the CESTAC method. With DSA, a real number becomes an N -dimensional set and any operation on these N -dimensional sets is performed element per element using the random rounding mode. The number of exact significant digits of such an N -dimensional set can be estimated from Eq. 1. The self-validation of the CESTAC method also leads to the concept of computational zero [START_REF]Zéro mathématique et zéro informatique[END_REF] defined below. Definition 2.1. During the run of a code using the CESTAC method, a result R is a computational zero, denoted by @.0, if ∀i,

R i = 0 or C R ≤ 0.
Any computed result R is a computational zero if either R = 0, R being significant, or R is insignificant. A computational zero is a value that cannot be differentiated from the mathematical zero because of its rounding error. From the concept of computational zero, an equality concept and order relations have been defined for DSA. Definition 2.2. Let X and Y be N -samples provided by the CESTAC method.

• Discrete stochastic equality denoted by ds= is defined as Xds= Y if and only if X -Y = @.0. Stochastic relational operators ensure that in a branching statement the same sequence of instructions is performed for all the samples which represent a variable. DSA enables to estimate the impact of rounding errors on any result of a scientific code and also to check that no anomaly occurred during the run, especially in branching statements.

Accuracy estimation by CADNA and SAM

The CADNA 1 software [START_REF] Jézéquel | CADNA: a library for estimating round-off error propagation[END_REF], [START_REF] Eberhart | High performance numerical validation using stochastic arithmetic[END_REF] is a library which implements DSA in any code written in C, C++ or Fortran and allows one to use new numerical types: the stochastic types. In essence, classical floating-point variables are replaced by the corresponding stochastic variables, which are composed of three perturbed floating-point values. The library contains the definition of all arithmetic operations and order relations for the stochastic types. The control of the accuracy is performed only on variables of stochastic type. When a stochastic variable is printed, only its exact significant digits appear. For a computational zero, the string "@.0" is printed. In contrast to interval arithmetic, that computes guaranteed results, the CADNA software provides, with the probability 95% the number of exact significant digits of any computed result.

The SAM library 2 [START_REF] Graillat | Stochastic Arithmetic in Multiprecision[END_REF] implements in arbitrary precision the features of DSA: the stochastic types, the concept of computational zero and the stochastic operators. The SAM library is written in C++ and is based on MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]. The particularity of SAM (compared to CADNA) is the arbitrary precision of stochastic variables. The SAM library with 24-bit (resp. 53-bit) mantissa length is similar to CADNA in single (resp. double) precision, except the range of the exponent is only limited by the machine memory. In SAM, the number of exact significant digits of any stochastic variable is estimated with the probability 95%, whatever its precision. Like in CADNA, the arithmetic and relational operators in SAM take into account rounding error propagation.

In CADNA and SAM all operators are overloaded, therefore their use in a program requires only a few modifications: essentially changes in the declarations of variables and in input/output statements. CADNA and SAM can detect numerical instabilities which occur during the execution of the code. Such instabilities are usually generated by numerical noise, i.e. a result having no correct digits.

3 Computation of multiple roots of polynomials using DSA

Newton method in stochastic arithmetic

In section 3.2, an algorithm to compute multiple roots of polynomials in stochastic arithmetic will be described. This algorithm uses Newton method whose behaviour in stochastic arithmetic has been studied in [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF]. We recall here properties of Newton method in stochastic arithmetic.

To compute the root α of a polynomial P , Newton method consists, from an initial approximation x 0 , in computing the following sequence for n ≥ 0:

x n+1 = x n - P (x n ) P (x n ) . ( 2
)
In classical floating-point arithmetic, Newton's iterations are stopped thanks to a criterion such as |x n+1 -x n | ≤ ε, but it may be difficult to choose a suitable value for the parameter ε taking into account rounding errors. In stochastic arithmetic, the stopping criterion becomes x n+1 ds= x n in accordance with Definition 2.2. It means iterations are stopped when the difference between two successive iterates is a computational zero. Thus useless iterations that would bring numerical noise to the solution are avoided. Algorithm 1 presents a stochastic version of Newton method that relies on this optimal stopping criterion.

Algorithm 1: stochastic Newton method: st-Newton Data: a polynomial P and an initial approximation x 0 of one of its roots Result: an approximation x of a root of P 1 x = x 0 ; 2 do 3 y = x; 4 x = y -P (y)/P (y); 5 while x ds = y;

The link between the exact root of a polynomial and the approximation provided by Algorithm 1 is discussed below. First, we need to make clear the notion of decimal significant digits in common between two numbers. Definition 3.1. The number of decimal significant digits in common between two real numbers a and b is defined in R by

• for a = b, C a,b = log 10 a + b 2(a -b) ; • for all a ∈ R, C a,a = +∞.
Then |a -b| = a+b 2 10 -C a,b . For instance, if C a,b = 3, the relative difference between a and b is of the order of 10 -3 , which means that a and b have three significant decimal digits in common.

The following theorem, introduced in [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF], gives a relation between the common significant digits of two successive approximations of the root computed using Newton method and the common significant digits of an approximation and the root itself, when the root is single. Theorem 3.1. Let x n and x n+1 be two successive approximations computed using Newton method of a polynomial root α of multiplicity m = 1 (single root).

Then

C xn,xn+1 ∼ ∞ C xn,α .
According to Theorem 3.1, in the convergence zone, the digits that are common to two successive approximations are also in common with the exact root. In stochastic arithmetic, Newton's iterations are stopped when x n+1 ds= x n , i.e. when x n+1 -x n is a computational zero (either numerical noise or the mathematical zero). In that case, the digits in x n+1 and x n that are not affected by rounding errors are the same and these digits are in common with the exact root. Therefore, according to Theorem 3.1, in the approximation provided by Algorithm 1, the digits not affected by rounding errors that are estimated thanks to DSA are those of the exact root.

Algorithm based on GCD and Newton method

When used for the computation of a single root of a polynomial, the convergence of Newton method is quadratic. But it becomes linear for the computation of a multiple root. In that case, the convergence speed of modified Newton method is satisfactory, but it requires the root multiplicity.

We propose Algorithm 2 that takes benefit of both the convergence speed of Newton method and stochastic arithmetic to control rounding errors. To compute the (possibly multiple) roots of a polynomial P , first we compute G, the stochastic GCD of the polynomial P and its derivative P : G = st-gcd(P, P ) using Algorithm 4. Then we compute the polynomial Q = P/G using a stochastic version of the Euclidean division (Algorithm 3). Algorithms 3 and 4 are stochastic versions of respectively the polynomial GCD and the polynomial Euclidean division. They use classic algorithms combined with the features of DSA, in particular the discrete stochastic relations recalled in Definition 2.2. The impact of DSA on these algorithms is described in Section 3.3.

The polynomial Q has normally only single roots and its degree d is its number of roots. Until a degree 4, its roots can be computed using adequate formulas, based on simple and well known expressions if d is 1 or 2, on Cardan's method if d is 3, and on Ferrari's method if d is 4 [START_REF] Kurosh | Higher algebra[END_REF]. So, if d ≤ 4, initial approximations of the roots are actually not required. From a degree 5, the stochastic version of Newton method (Algorithm 1) is used. The polynomial Q having single roots, Newton method exhibits a quadratic convergence.

An advantage of Algorithm 2 is the fact that the polynomial Q is computed once whatever the number of roots. Then all the roots are computed, either thanks to adequate formulas, or using Newton method applied for each root to the same low-degree polynomial.

Benefits of the DSA implementation

We exemplify the benefits of DSA with the computation of the root of a low-degree polynomial: P (x) = (3x -1) 5 . The first step in Algorithm 2 is the computation of G = st-gcd(P, P ) using Algorithm 4. The polynomial R 2 in Algorithm 4 should be null, and the returned result should be R 1 = P . Then, the polynomials G and Q in Algorithm 2 should be respectively of degree 4 and 1.

We compare two kinds of executions:

• with Algorithms 2, 3 and 4 using DSA and in particular stochastic relations recalled in Definition 2.2;

Algorithm 2: Computation of polynomial roots based on st-gcd and st-Newton Data: a polynomial P and an array X 0 of initial approximations of its roots Result: an array X of approximations of the roots 1 G = st-gcd(P, P ); 2 Q = quotient(st-Euclidean-div(P , G)); 3 if degree(Q) ≤ 4 then 4 computation of X using adequate formulas; 5 end 6 else

7 for i=1 to degree(Q) do 8 X[i]=st-Newton(Q, X 0 [i]);
≤ m ≤ n Result: polynomials Q and R s.t. A = B * Q + R 1 R = A; 2 for i = n -m to 0 do 3 if lc(R) ds = 0 then // tests if degree(R) = m + i // lc(R): leading coefficient of R 4 q i = lc(R)/lc(B); 5 R = R -q i x i B;
1 R 0 = A; 2 R 1 = B; 3 i = 1; 4 while R i ds = 0 do 5 R i+1 =remainder (st-Euclidean-div(R i-1 , R i )); // Ri+1 = Ri-1 Mod Ri 6 i = i + 1; 7 end 8 return R i-1 ;
• with the same algorithms implemented using classic floating-point arithmetic. Because we focus on partial results obtained at the first steps of the computation, naive comparisons to zero are performed in Algorithms 3 and 4.

Table 1 presents results computed with/without DSA using a working precision of 35, 36 or 37 bits. We report the coefficients of the polynomial R 2 (x) = 3 i=0 R i 2 x i , the degree of the polynomials G and Q, and the computed root.

Without DSA, the computation is carried out using MPFR. Because we are particularly interested in the first iteration in Algorithm 4, the condition in line 4 is a naive comparison to zero. With 35 or 37 bits, all the coefficients of R 2 are zeroes, therefore the degrees of polynomials G and Q are correct. An approximation of the root is computed, but we have no information on its numerical quality. With 36 bits, three coefficients in R 2 are non-zeroes because of rounding errors. Undesirable iterations are performed in Algorithm 4 and consequently polynomials G and Q have incorrect degrees. It may be difficult to find a suitable stopping criterion in GCD computation that discards coefficients affected by rounding errors, but not small coefficients correctly computed. Furthermore the coefficients in R i may be numerical noise with different orders of magnitude. For instance if the root of (3x-1) 10 is computed using Algorithm 2 with 36 bits, because of rounding errors, R 2 has non-zero coefficients with orders of magnitude that vary from 10 -6 to 10 -10 .

With the DSA implementation based on the SAM library, the condition in Algorithm 4 (on line 4) is satisfied if all the coefficients of the polynomial R i are different from a computational zero in accordance with Definitions 2.1 and 2.2. This stopping criterion does not require any ε parameter and enables one to discard numerical noise whatever the working precision. Furthermore DSA estimates which digits in the results are not affected by rounding errors and displays only these correct digits. One can observe in Table 1 that the coefficients of the polynomial R 2 are identified as numerical noise, the degrees of G and Q are correctly computed, and the root is provided with 9 digits estimated as correct by DSA.

Algorithms 1 to 4 benefit from various features of DSA.

•

In Algorithm 1, line 5, as already mentioned in Section 3.1, the criterion enables one to stop when the difference between two successive iterates is a computational zero, and consequently useless iterations are avoided.

•

In Algorithm 2, if the degree d of the polynomial Q satisfies d ≤ 4, its roots are directly computed using arithmetic expressions and the rounding errors affecting them are estimated by DSA. If d > 5, because Q has single roots, as already mentioned in Section 3.1, in each computed root, the digits estimated by DSA as not affected by rounding errors are those of the exact root.

•

In Algorithm 3, line 3, the test is achieved by checking whether the leading coefficient of R is a computational zero.

0 0 -1.86264514923e-9 0 @.0 R 2 2 0 0 9.31322574615e-10 0 @.0 R 1 2 0 0 -2.32830643654e-10 0 @.0 R 0 2 0 0 0 0 @.0 degree(G) 4 4 0 4 4 degree(Q) 1 1 5 1 1 root 1/3 3.33333333328e-1 Fail 3.333333333321e-1 0.333333333 TABLE 1
Computation of the root of (3x -1) 5 using Algorithm 2 with/without DSA

Numerical experiments

Numerical experiments have been carried to cover two aspects: the accuracy of the computed results and the performance. The platform for performance measurements is an Intel Core i7-8650U processor clocked at 1.9 GHz with 8 MB cache. The codes are compiled with gcc version 8.4.0 and optimized with the -O3 flag. The instability detection level chosen in the SAM library is the self-validation, already introduced in 2.1: all multiplications and divisions are controlled during the execution.

The roots of various polynomials have been computed using Algorithm 2 with the SAM library. The results have been compared with those also computed with SAM, but using an algorithm described in [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] and summarized below. In [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] a root of a polynomial P with multiplicity m > 1 is evaluated using modified Newton method that consists, from an initial approximation x 0 , in computing the following sequence for n ≥ 0:

x n+1 = x n -m P (x n ) P (x n ) . ( 3
)
For multiple polynomial roots, modified Newton method exhibits a better convergence than Newton method. However it requires the evaluation of the multiplicity m. In [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] the multiplicity is determined thanks to three successive iterations of Newton method as proposed in [START_REF] Yakoubsohn | Numerical elimination, Newton method and multiple roots[END_REF]. Then the root is approximated using modified Newton method. In Algorithm 2 the working precision depends on the requested accuracy and a rate both chosen by the user:

P recision = Requested_accuracy * Rate with Rate > 1. (4)
Indeed because of rounding errors, the working precision has to be greater than the requested accuracy. In [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] the initial precision is set according to Eq. 4. Then the working precision is doubled between successive calls to modified Newton method applied iteratively.

Numerical experiments have been carried out with two types of polynomials.

Q n = (3x-2) n1 (7x-3) n2 (13x-4) n3 (19x-2) n4 (23x- 1) n5
The roots of Q n are denoted as

β 1 = 2/3, β 2 = 3/7, β 3 = 4/
n = 4 i=1 n i = 5000 with n i = 1000 (i = 1, • • • , 5).
Table 2 presents results obtained using Algorithm 2, Tables 3 and 4 results computed using modified Newton method. As a recall, Algorithm 2 computes all the roots of a polynomial, whereas the algorithm based on modified Newton method requires an execution per root. Table 2 presents for polynomials P 54 , P 104 , Q 55 , Q 105 , and Q 5000 :

• the requested number of decimal digits in all the roots,

• the number of digits in common with the exact roots; this number depending on the root, the minimum and the maximum number of digits are mentioned,

• the minimum rate that enables the accuracy requirement to be fulfilled for all the roots; this rate is determined starting from 1.1 and increasing by steps of 0.1,

• the execution time in seconds,

• the performance ratio w.r.t. the algorithm from [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] based on modified Newton method. In Algorithm 2 after the execution of st-gcd and st-Euclidean-div, the four roots of polynomials P 54 and P 104 are computed using Ferrari's method [START_REF] Kurosh | Higher algebra[END_REF], whereas the five roots of Q 55 , Q 105 , and Q 5000 are computed using st-Newton. It can be observed in Table 2 that the accuracy requirement is fulfilled with a multiplicative rate between 1.1 and 1.5. For each root, the number of digits in common with the exact result is also the number of digits not affected by rounding errors estimated by DSA. This equality is in accordance with Theorem 3.1. Algorithm 2 outperforms the algorithm based on modified Newton method by 1 to 4 orders of magnitude. Concerning the roots of P 104 with 1000 digits, the performance ratio is not available because the algorithm based on Newton method has failed to compute one of the roots. It must be pointed out that the first 1000 digits of the roots of polynomial Q 5000 (of degree 5000) are computed in less than 1 second and the execution time becomes 5 seconds if 5000 digits are required.

Tables 3 and4 present for each root of the polynomials P 54 , P 104 , Q 55 , and Q 105 :

• the requested number of decimal digits in the root,

• the number of digits not affected by rounding errors estimated by DSA,

• the number of digits in common with the exact root,

• the minimum rate that enables the accuracy requirement to be fulfilled; again this rate is determined starting from 1.1 and increasing by steps of 0. It can be observed in Tables 3 and4 that, for each root, the number of digits not affected by rounding errors estimated by DSA is greater than the number of digits in common with the exact result. The ratio between these numbers varies from 1.3 to 12. As a remark, Theorem 3.1 applies to the computation of roots of multiplicity 1 using Newton method, not to the computation of multiples roots using modified Newton method. The minimum multiplicative rate for the accuracy requirement to be fulfilled varies from 1.6 to 9.7 depending on the requested accuracy and the root multiplicity.

In particular, a rate of 9.7 is required to obtain the first 500 digits of the root α 4 of P 104 that has a multiplicity of 50. Modified Newton method fails to compute the same root with 1000 digits. With a rate of 10.9, the method provides the first 984 digits of α 4 . But the working precision associated to a rate of 11, that is initially 11,000 digits and then increases, results in a failure due to the memory limit. For all root approximations based on modified Newton method, the root multiplicity is correctly determined. However, then modified Newton method is used with a possibly high-degree polynomial (from degree 54 to 105 in our experiments). With Algorithm 2, Newton method is used with a low-degree polynomial (a polynomial of degree 4 or 5 for the results presented in Table 2).

Conclusion and perspectives

In this paper, we have presented an algorithm that efficiently and accurately computes polynomial roots, in particular multiple roots. This algorithm uses stochastic arithmetic that enables one to estimate rounding errors, and so to discard results that are actually numerical noise. With stochastic arithmetic, iterative algorithms can be stopped in an optimal way that does not rely on any parameter. Thanks to a stochastic version of the polynomial GCD algorithm and the polynomial Euclidean division, the proposed algorithm provides a low-degree polynomial with single roots. Depending on the number of roots, they can be computed using either adequate formulas or Newton method. In the computed roots, stochastic arithmetic estimates which digits are not affected by rounding errors, and these digits are in common with the exact roots. W.r.t. an algorithm presented in [START_REF] Graillat | Dynamical Control of Newton's Method for Multiple Roots of Polynomials[END_REF] and based on modified Newton method, the proposed algorithm requires less working precision and less execution time: in our experiments the performance gain is up to 4 orders of magnitude.

The proposed algorithm takes into account an accuracy requirement and a multiplicative rate to determine the required working precision. From our experiments, this rate remains low and varies from 1.1 to 1.5. However it would be interesting to determine it automatically from information such as the required accuracy, the polynomial degree, and the multiplicity of the roots. Another perspective would be, for polynomials having at least 5 roots, the computation in parallel of their roots. Indeed, in this case, the roots are computed thanks to independent executions of Newton method.

3 :

 3 Stochastic polynomial Euclidean division: st-Euclidean-divData: polynomial A of degree n, polynomial B of degree m, with 0

Algorithm 4 :

 4 q i x i and R stochastic polynomial GCD: st-gcd Data: polynomials A and B Result: stochastic GCD of A and B

•

  P n = (19x + 5) n1 (19x + 21) n2 (19x + 46) n3 (19x + 67) n4 The roots of P n are denoted as α 1 = -5/19, α 2 = -21/19, α 3 = -46/19, and α 4 = -67/19. The degrees of polynomials P n have been chosen as follows. n = 4 i=1 n i = 54 with n 1 = 7, n 2 = 9, n 3 = 13, n 4 = 25 n = 4 i=1 n i = 104 with n 1 = 10, n 2 = 18, n 3 = 26, n 4 = 50.

•

5 i=1

 5 13, β 4 = 2/19, and β 5 = 1/23. The degrees of polynomials Q n have been chosen as follows. n = n i = 55 with n 1 = 13, n 2 = 12, n 3 = 11, n 4 = 10, n 5 = 9 n = 4 i=1 n i = 105 with n 1 = 18, n 2 = 19, n 3 = 21, n 4 = 22, n 5 = 25

Acknowledgment

This work was supported by the InterFLOP (ANR-20-CE46-0009) project of the French National Agency for Research (ANR).