Improving the Performance of the Continued Fractions Method Using New Bounds of Positive Roots

Abstract. In this paper we compare four implementations of the Vincent -AkritasStrzeboński Continued Fractions ( VAS-CF) real root isolation method using four different (two linear and two quadratic complexity) bounds on the valu es of the positive roots of polynomials. The quadratic complexity bounds were inclu ded to see if the quality of their estimates compensates for their quadratic complex ity. Indeed, experimentation on various classes of special and random polynomials reveal d that theVAS-CF implementation usingLMQ, theQuadratic complexity variant of our LocalMax bound, achieved an overall average speed-up of 40 % over the original implementation using Cauchy’s linear bound.


Introduction
We begin by first reviewing some basic facts about the VAS-CF method for isolating the positive roots of polynomials.This method is based on Vincent's theorem of 1836, [1], which states: Theorem 1.If in a polynomial, p(x), of degree n, with rational coefficients and without multiple roots we perform sequentially replacements of the form where α 1 ≥ 0 is an arbitrary non negative integer and α 2 , α 3 , . . .are arbitrary positive integers, α i > 0, i > 1, then the resulting polynomial either has no sign variations or it has one sign variation.In the last case the equation has exactly one positive root, which is represented by the continued fraction , whereas in the first case there are no positive roots.
See the papers by Alesina & Galuzzi, [2] and Chapter 7 in [3] for a complete historical survey of the subject and implementation details respectively.The thing to note is that the quantities α i (the partial quotients of the continued fraction) are computed by repeated application of a method for estimating lower bounds 1 on the values of the positive roots of a polynomial.
Cauchy's (linear complexity) bound on the values of the positive roots of a polynomial, was used until recently in the VAS-CF real root isolation method, [4].In the SYNAPS implementation of the VAS-CF method, [5], Emiris and Tsigaridas used Kioustelidis' (linear complexity) bound, [6] and independently verified the results obtained by Akritas and Strzeboński, [4] 2 .
In this paper we implemented the VAS-CF real root isolation process using both linear and quadratic complexity bounds on the values of the positive roots of a polynomial -all of which are based on Theorem 2 below, [9], [10], [11] -and discovered that the latter substantially improve its performance!
The rest of the paper is structured as follows: In Section 2 we present the VAS-CF algorithm, as found in [4]; this is done both for completeness of the present paper and to correct a misprint that appeared in Step 5 in [4].We also present the theoretical background of our new quadratic complexity bounds.
In Section 3, we compare four different implementations of the VAS-CF algorithm using two linear and two quadratic complexity bounds -all suitably adjusted for computing lower bounds on the values of the positive roots.
Finally, in Section 4 we present our conclusions; namely, when LMQ, the Quadratic complexity variant of our Local Max bound, is implemented in VAS-CF, an overall average speed-up of 40 % is achieved over the original version of VAS-CF, where Cauchy's linear bound is implemented.

Algorithmic and theoretical background
In Subsection 2.1 we present the VAS-CF algorithm -as found in [4] -and correct a misprint in Step 5 that had appeared in that presentation; moreover, we explain where the new bound on the positive roots is used.We also present our Theorem 2 from which all the bounds on the values of the positive roots are derived.Linear complexity bounds are presented in Subsection 2.2, whereas the quadratic complexity ones are presented in Subsection 2.3.

Description of the Continued Fractions Algorithm VAS-CF
Using the notation of paper [4], let f ∈ Z[x] \ {0}.By sgc(f ) we denote the number of sign changes in the sequence of nonzero coefficients of f .For nonnegative integers a, b, c, and d, such that ad − bc = 0, we put The value of parameter α 0 used in Step 4 below needs to be chosen empirically.In our implementation α 0 = 16.
3. Compute a lower bound α ∈ Z on the positive roots of p.
Please note that the lower bound, α, on the positive roots of p(x) is computed in Step 3, and used in Step 5.
As mentioned in the introduction, Cauchy's bound was the first one to be used in VAS-CF.However, new bounds on the values of the positive roots were developed by us recently; for details see [9], [10], and [11].The important thing that emerged from our work is that all bounds on the values of the positive roots are derived from the following [11]: be a polynomial with real coefficients and let d(p) and t(p) denote the degree and the number of its terms, respectively.Moreover, assume that p(x) can be written as where all the polynomials q i (x), i = 1, 2, . . ., 2m and g(x) have only positive coefficients.In addition, assume that for i = 1, 2, . . ., m we have and , where e 2i−1,1 = d(q 2i−1 ) and e 2i,1 = d(q 2i ) and the exponent of each term in q 2i−1 (x) is greater than the exponent of each term in q 2i (x).If for all indices i = 1, 2, . . ., m, we have then an upper bound of the values of the positive roots of p(x) is given by for any permutation of the positive coefficients c 2i−1,j , j = 1, 2, . . ., t(q 2i−1 ).Otherwise, for each of the indices i for which we have t(q 2i−1 ) < t(q 2i ), we break up one of the coefficients of q 2i−1 (x) into t(q 2i ) − t(q 2i−1 ) + 1 parts, so that now t(q 2i ) = t(q 2i−1 ) and apply the same formula (3) given above.
For a proof of this theorem see [11].Among others, the following linear and quadratic complexity bounds on the values of the positive roots of polynomials are derived from it.

Linear complexity bounds derived from Theorem 2
Various linear complexity bounds can be obtained from the above theorem; the ones described below have been presented elsewhere, [11], but not in the context of complexity.We present them here again, briefly, for completeness: C. Cauchy's "leading-coefficient" implementation of Theorem 2. For a polynomial p(x), as in equation ( 1), with λ negative coefficients, Cauchy's method first breaks up its leading coefficient, α n , into λ equal parts and then pairs each part with the first unmatched negative coefficient, [12].
That is, we have: or, equivalently, K. Kioustelidis' "leading-coefficient" implementation of Theorem 2. For a polynomial p(x), as in (1), Kioustelidis' method matches the coefficient −α n−k of the term −α n−k x n−k in p(x) with αn 2 k , the leading coefficient divided by 2 k .Kioustelidis' "leading-coefficient" implementation of Theorem 2 differs from that of Cauchy's only in that the leading coefficient is now broken up in unequal parts, by dividing it with different powers of 2, [6], or, equivalently, .
In addition to the above we present our own two linear complexity implementations of Theorem 2, the combination of which (that is, taking the minimum) yields the best linear complexity upper bound on the values of the positive roots of a polynomial, [11].Please note that in general we can obtain better bounds if in Theorem 2 we pair coefficients from non adjacent polynomials q 2ℓ−1 (x) and q 2i (x), 1 ≤ ℓ < i -and this is done in the following definitions: FL. "First-λ" implementation of Theorem 2. For a polynomial p(x), as in (2), with λ negative coefficients we first take care of all cases for which t(q 2i ) > t(q 2i−1 ), by breaking up the last coefficient c 2i−1,t(q2i) , of q 2i−1 (x), into t(q 2i ) − t(q 2i−1 ) + 1 equal parts.We then pair each of the first λ positive coefficients of p(x), encountered as we move in non-increasing order of exponents, with the first unmatched negative coefficient.
LM. "Local-Max" implementation of Theorem 2. For a polynomial p(x), as in (1), the coefficient −α k of the term −α k x k in p(x) -as given in (1) -is paired with the coefficient αm 2 t , of the term α m x m , where α m is the largest positive coefficient with n ≥ m > k and t indicates the number of times the coefficient α m has been used.
We have tested extensively -on various classes of specific and random polynomialsall four linear complexity bounds mentioned above and the following is a summary of our findings, [11]: • Kioustelidis' bound is, in general, better (or much better) than Cauchy's; this happens because the former breaks up the leading coefficient in unequal parts, whereas the latter breaks it up in equal parts.
• Our first-λ bound, as the name indicates, uses additional coefficients and, therefore, it is not surprising that it is, in general, better (or much better) than both previous bounds.In the few cases where Kioustelidis' bound is better than first-λ, our localmax bound takes again the lead.
Therefore, given their linear cost of execution, min(FL, LM) is the best among the linear complexity bounds, [11], and we use it in Section 3 to determine the cutoff point between linear and quadratic complexity bounds.Additionally, in Section 3 we see that implementing min(FL, LM) in the continued fractions real root isolation method, VAS-CF/ min(FL, LM), we obtain a speed-up of 15 % -when compared with V AS − CF/Cauchy.

Quadratic complexity bounds derived from Theorem 2
In this subsection we present the two quadratic complexity bounds which will be used in our experimentation.See the literature, [13], for a more complete discussion of quadratic complexity bounds.We begin with KQ, the quadratic version of Kioustelidis' bound on the values of the positive roots of a polynomial.To our knowledge this was first presented by Hong, [14, p. 574], and can be stated as follows: KQ.Kioustelidis' Quadratic complexity implementation of Theorem 2. For a polynomial p(x), as in equation ( 1), each negative coefficient a i < 0 is "paired" with each one of the preceding positive coefficients a j divided by 2 j−i -that is, each positive coefficient a j is "broken up" into unequal parts, as is done with just the leading coefficient in Kioustelidis' bound -and the minimum is taken over all j; subsequently, the maximum is taken over all i.
That is, we have: or, equivalently, Hong's Theorem 2.2, if examined carefully, reveals that for polynomials in one variable it is the quadratic version of Kioustelidis' bound, KQ; moreover, the relation with our Theorem 2 is quite obvious.This observation lead us to our own quadratic complexity bound, LMQ, based also on our Theorem 2, [13], which can be stated as follows: LMQ. "Local-Max" Quadratic complexity implementation of Theorem 2. For a polynomial p(x), as in (1), each negative coefficient a i < 0 is "paired" with each one of the preceding positive coefficients a j divided by 2 tj -that is, each positive coefficient a j is "broken up" into unequal parts, as is done with just the locally maximum coefficient in the local max bound; t j is initially set to 1 and is incremented each time the positive coefficient a j is used -and the minimum is taken over all j; subsequently, the maximum is taken over all i That is, we have: Since 2 tj ≤ 2 j−i -where i and j are the indices realizing the max of min; equality holds when there are no missing terms in the polynomial -it is clear that the bounds computed by LMQ are sharper by the factor 2 j−i−t j j−i than those computed by Kioustelidis' KQ; nonetheless, VAS-CF is implemented with both quadratic bounds.
We have also implemented the fastest quadratic complexity bound, F LQ, but we omit it since its behavior is similar to that of LMQ and it cannot be directly compared with KQ; details can be found elsewhere, [15].

Empirical results
In this section we compare four implementations of the VAS-CF real root isolation method using two linear and two quadratic complexity bounds on the values of the positive roots of polynomials.
The two linear complexity bounds are: Cauchy's and min(FL, LM), the minimum of our First Lambda and Local Max bounds, [11], whereas the two quadratic complexity ones are: KQ, the Quadratic complexity variant of Kioustelidis' bound, studied by Hong, [14], and LMQ, the Quadratic complexity version of our Local Max bound, [13].
Our choice of the various bounds in the implementations of VAS-CF is justified as follows: 1. From the linear complexity bounds we included: (a) Cauchy's bound, C, to be used as a point of reference, since it has been in use for the past 30 years, and (b) our min(FL, LM) bound, [11], which is the best among the linear complexity bounds, in order to see when it's implementation will outperform that of the two quadratic complexity bounds.
2. From the quadratic complexity bounds we included: (a) Kioustelidis' KQ, and (b) our own LMQ, in order to compare their performance; as explained in the previous section LMQ computes sharper estimates than KQ.
We have implemented all versions of VAS-CF as a part of Mathematica kernel.They all use the same implementation of Shaw and Traub's algorithm for Taylor shifts (see [16]).We followed the standard practice and used as benchmark the Laguerre 3 , Chebyshev (first 4 and second 5 kind), Wilkinson 6 and Mignotte 7 polynomials, as well as several types of randomly generated polynomials of degrees {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1500, 2000}.For the random polynomials the size of the coefficients ranges from −2 20 to 2 20 .
All computations were done on a Windows XP laptop computer with with 1.8 Ghz Pentium M processor, and 2 GB of RAM.The change in memory use was negligible in every case, and hence, it is not included in the following tables.The average speed-up was calculated using the formula: Speed-up = 100 • (LMQ − Cauchy)/Cauchy, from which we omitted the minus sign. 3recursively defined as: 4 recursively defined as: T 0 (x) = 1, T 1 (x) = x, and T n+1 (x) = 2xTn(x) − T n−1 (x) 5 recursively defined as: by interval data we denote a list {a, b, c, d, p, s}, where p is a polynomial such that the roots of f in intrv(a, b, c, d) are images of positive roots of p through Φ a,b,c,d , and s = sgc(p).

Table 1 .
Special Polynomials.Some indicative results are given for degrees 100, 1000, 1500 and 2000.The average speed-up for this table is about 32 %

Table 3 .
Polynomials with random 1000-bit coefficients.The average speed-up for this table is about 48 %

Table 4 .
Monic polynomials with random 10-bit coefficients.The average speed-up for this table is about 31 %

Table 5 .
Monic polynomials with random 1000-bit coefficients.The average speed-up for this table is about 60 %