Effective Real Root Isolation Using Continued Fractions

Alkiviadis G. Akritas

Adam W. Strzebonski

Recent progress in polynomial elimination has rendered the computation of the real roots of ill-conditioned polynomials of high degree > 1000 with huge coefficients (several thousand digits) a critical operation in computer algebra. To rise to the occasion, the only method-candidate that has been considered by various authors for modification and improvement has been the Collins-Akritas bisection method [CA1976]. The most recent example is a preprint [RZ2002], where the authors try to present

...all the known algorithms in a unified framework. Using that framework, a new algorithm is presented, which is optimal in terms of memory usage and as fast as both Collins and Akritas' algorithm and Krandick variant... 

However, the authors of that preprint failed to include our own continued fractions method [ABS1994]  in their consideration. Developed in 1994 for implementation in the computer algebra system Mathematica, our method is based on Vincent's theorem [V1836]-extended by Uspensky [U1948], to guarantee that the process terminates, and modified by one of us [Akr1978]-and is amply described in [Akr1980], [Akr1989] and [ABS1994]. Experimentation with the data presented in [RZ2002] showed that, with respect to time, our continued fractions method is by far superior to the one presented in the most recent preprint, whereas the two are about equal with respect to space. Moreover, it should be kept in mind that there is only one unifying principle, to wit continued fractions, which in the Collins-Akritas bisection method are indirectly computed from the intervals. Therefore, our method is the unifying framework.

Introduction-or Paying tribute to Vincent!

Vincent's theorem, extended to guarantee termination, can be stated as follows:

Vincent's Theorem (1836). Let p(x) = 0 be a polynomial equation of degree n > 1, with rational coefficients and without multiple roots, and let Δ > 0 be the smallest distance between any two of its roots. Let m be the smallest index such that

(Subscript[F, m - 1]) (Δ/2) > 
      1, (Subscript[F, m - 1]) (Subscript[F, m]) Δ > 
      1 + (1/(Subscript[ϵ, n])),

where (Subscript[F, k]) is the k-th member of the Fibonacci sequence 1, 1, 2, 3, 5, 8, 13, 21, … and

(Subscript[ϵ, n]) = (((1 + (1/n)))^((1/(n - 1)))) - 
        1

Let (Subscript[a, 1]) be an arbitrary nonnegative integer, and let  (Subscript[a, 2]), …, (Subscript[a, m]) be arbitrary positive integers. Then the substitution

x⟵(Subscript[a, 1]) + 
      (1/((Subscript[a, 2]) + (1/(⋱ + (1/((Subscript[a, m]) + (1/ξ)))))))

transforms the equation  p(x) = 0  into the equation p(ξ) = 0, which has not more than one sign variation.

This theorem can be used to isolate the positive roots of a polynomial equations p(x) = 0. For the negative roots we simply replace x⟵-x and repeat the process. The substitution of x, in p(x) = 0, by the continued fraction shown in Vincent's theorem results in all roots being clustered around -1 within a circle of radius (Subscript[ϵ, n]).

Vincent used the weaker version of his theorem (that is, without the terminating condition) to isolate the real roots of polynomial equations with rational coefficients [V1836]. In his paper of 1836 he presents several examples to demonstrate the concepts involved.

Note that Vincent's method appeared several years after Sturm's method for isolating the real roots. Due to Sturm's fame and his priority, Vincent's method was almost totally forgotten.

In 1975-76 Akritas rediscovered Vincent's theorem in Uspensky's book. Without thoroughly studying the historical background, Collins and Akritas-in their paper [CA1976]-erroneously named the method they developed, their method, the modified Uspensky's method. They were misled by Uspensky's claim (in the preface of his book [U1948]) that he himself had invented Vincent's method.

The mistake was caught by Akritas in his Ph.D.Thesis [Akr1978b] and was subsequently announced in [Akr1982], [Akr1984] and [Akr1989]. The gist of the matter is that Uspensky was not aware of the Budan's theorem; its statement is as follows:

Budan's Theorem (1807). If in an equation in x of degree n > 0, p(x) = 0 we make two substitutions, x⟵α + x' and x⟵ β + x'', where α and β are real numbers such that α < β, then the following hold:

(α) The transformed equation in x' = x - α cannot have fewer sign variations than the transformed equation in x'' = x - β.

(β) The number of real roots in of the equation p(x) = 0, located between α and β, can never be greater than the number of sign variations lost in the sequence of coefficients in passing from the transformed equation in x' = x - α to the transformed equation in x'' = x - β.

(γ) When the number of real roots of p(x) = 0 located between α and β is less than the number of sign variations lost in the sequence of coefficients in passing from the transformed equation in x' = x - α to the transformed equation in x'' = x - β, the difference is an even number.

Note that Budan's theorem above is almost always confused with the one by Fourier on the same subject. It is exclusively Fourier's theorem of 1820 that appears in the literature under the names "Fourier", "Budan", or "Budan-Fourier".

Because of his ignorance of Budan's theorem, Uspensky uses Vincent's theorem in the following time consuming way to compute a partial quotient (Subscript[a, i]): after each translation (Taylor shift) x⟵x + 1 he has to perform the substitution x⟵(1/(x + 1)) , to make sure he has not lost any roots in that interval.

Vincent, on the contrary, being fully aware of Budan's theorem-he actually states both Budan's and Fourier's theorems in his 1836 article-makes only Taylor shifts x⟵x + 1 until he detects a loss in the number of sign variations. Then and only then does Vincent make the substitution x⟵(1/(x + 1)) knowing that he has computed a partial quotient (Subscript[a, i]). Therefore, what Uspensky did was to take Vincent's method and double its computing time.

With this in mind it is somewhat baffling that the authors in [RZ2002] time and again refer to a certain Uspensky's method and keep referring to the Collins-Akritas method as the modified Uspensky's method. Hopefully this practice will end. It sets a very bad precedent!

The only problem with Vincent's continued fractions method is that it is exponential in the following two ways:

(exp\(TraditionalForm\`\_1\))  in the number of substitutions of the form x⟵x + 1, which need to be performed in order to compute
a partial quotient (Subscript[a, i]), and

(exp\(TraditionalForm\`\_2\)) in the growth of the integer coefficients of the polynomials, which are obtained after substitutions of the form x⟵x + (Subscript[a, i]).

One can easily create examples where the number of translations makes it prohibitive to use. Uspensky realized the problem but failed to solve it.

Our continued fractions method

It was in Akritas' Ph.D.Thesis [Akr1978b] that the exponentiality problem of Vincent's continued fractions method was solved; see also [Akr1980] and [Akr1989].

To tackle (exp\(TraditionalForm\`\_1\)), the problem with the number of transformations, Akritas used Cauchy's rule for computing a lower bound on the values of the positive roots of a polynomial equation. That way, instead of performing transformations of the form x⟵x + 1, he performs Taylor shifts of the form x⟵x + (Subscript[a, i]), where (Subscript[a, i]) = floor(nextroots). It should be noted that the substitutions  x⟵x + 1 and x⟵x + (Subscript[a, i]), for (Subscript[a, i]) > 1, are executed in about the same amount of time [AD1980].

To tackle (exp\(TraditionalForm\`\_2\)), the problem with the growth of the integer coefficients, Akritas used a theorem by Lang and Trotter [LT1972] to justify the claim that \(TraditionalForm\`a\_i = O(\(||\)\(p(x)\)\( || \_∞\))\) for all  i = 1, …, m; that is, the (Subscript[a, i])s are about as big as the biggest coefficient of the polynomial equation whose roots we want to isolate.

The interesting result by Lang and Trotter states that:

In the continued fraction expansion of almost all numbers, the probability that (Subscript[a, i]), the i-th partial quotient, is equal to a  positive integer j is given by (Subscript[log, 2]) ((((j + 1))^2)/(j(j + 2))).

For almost all numbers, this means that the probability for (Subscript[a, i]) = 1 is approximately 0.41. Likewise, the probability for (Subscript[a, i]) = 1000 is approximately 1.44 (10^(-6)) ≈ 0.

So, the result by Lang and Trotter tells us, on one hand, that we cannot bound the size of the partial quotients and, hence, we cannot control the growth of the coefficients. On the other hand, however, due to the fact that we use few partial quotients to isolate the roots, this same result tells us that the probability of big partial quotients is equal to zero and, hence, the claim is well-founded [Akr1989].

Subsequently Akritas, Bocharov and Strzebonski improved the performance of continued fraction method for polynomials with large roots by adding a substitution x⟵a x whenever the lower bound a on the values of the positive roots is large. In [ABS1994] they run extensive empirical comparison between implementations of the
continued fractions method and the Collins-Akritas method, and show that the continued fractions method performs significantly faster in all cases except for polynomials with all roots real and extermely large, where it is slightly slower.

Our recent improvement is to reduce the memory usage by the continued fraction method. We show that the algorithm can be so structured that the maximal number of transformed versions of the polynomial that need to be stored at any given time is bounded by 1 + (Subscript[log, 2]) n, where n is the degree of the input polynomial. This bound is based on the following conjecture, which we have not proven, but which we have extensively tested.

Let \(TraditionalForm\`f ∈ ℤ[x] \0\) be a polynomial of degree n, and let sgc(f) denote the number of sign changes in the sequence of nonzero coefficients of f.

Conjecture.

sgc (f (x + 1)) + 
        sgc ((((x + 1))^n) f ((1/(x + 1)))) ≤ 
      sgc (f (x))

Description of the algorithm

Let us first introduce the notation used in the algorithm. Let \(TraditionalForm\`f ∈ ℤ[x] \0\). For nonnegative integers a,b, c, and d such that a d - b c ≠ 0, we put

intrv (a, b, c, 
          d) := (Subscript[Φ, a, b, c, 
            d]) ((0, ∞))

where

(Subscript[Φ, a, b, c, 
          d]) : (0, ∞) ∋ 
        x⟶((a x + b)/(c x + 
                d)) ∈ (0, ∞)

and 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 (Subscript[Φ, a, b, c, d]), and s = sgc(p).

The value of parameter (Subscript[α, 0]) used in step 4 below needs to be chosen empirically. In our implementation (Subscript[α, 0]) = 16.

Algortihm Continued Fractions.
Input: A squarefree polynomial \(TraditionalForm\`f ∈ ℤ[x] \0\).
Output: The list  rootlist of positive roots of f.

(1) Set rootlist to an empty list. Compute s ← sgc (f).  If s = 0 return an empty list. If s = 1 return (0,∞). Put interval data {1, 0, 0, 1, f, s} on intervalstack.

(2) If  intervalstack is empty, return rootlist, else take interval data {a, b, c, d, p, s} off intervalstack.

(3) Compute a lower bound α on positive roots of p.

(4) If  α > (Subscript[α, 0]) set p(x)⟵p(α x), a ← α a, c ← α c,  α ← 1.

(5) If α ≥ 1, set  p(x)⟵p(x + α), b ← α a + b,and c ← α c. If p( 0) = 0, add [b/d, b/d] to rootlist, and set p(x)⟵p(x)/x. Compute s ← sgc (p). If s = 0 go to step 2. If s = 1, add intrv(a, b, c, d) to rootlist and go to step 2.

(6) Compute (Subscript[p, 1])(x)⟵p(x + 1), and set  (Subscript[a, 1]) ← a, (Subscript[b, 1]) ← a + b, (Subscript[c, 1]) ← c, (Subscript[d, 1]) ← c + d, and r ← 0. If (Subscript[p, 1])( 0) = 0, add [((Subscript[b, 1])/(Subscript[d, 1])), ((Subscript[b, 1])/(Subscript[d, 1]))] to rootlist, and set (Subscript[p, 1])(x)⟵(Subscript[p, 1])(x)/
        x, and r ← 1. Compute  (Subscript[s, 1]) ← sgc ((Subscript[p, 1])), and set  (Subscript[s, 2]) ← s - (Subscript[s, 1]) - r,  (Subscript[a, 2]) ← b, (Subscript[b, 2]) ← a + b, (Subscript[c, 2]) ← d, (Subscript[d, 2]) ← c + d.

(7) If (Subscript[s, 2]) > 1, compute (Subscript[p, 2])(x)⟵(((x + 1))^m)  p(
          (1/(x + 1))), where m is the degree of p. If (Subscript[p, 2])( 0) = 0, set (Subscript[p, 2])(x)⟵(Subscript[p, 2])(x)/
        x. Compute (Subscript[s, 2]) ← sgc ((Subscript[p, 2])).

(8) If (Subscript[s, 1]) < (Subscript[s, 2]), swap {(Subscript[a, 1]), (Subscript[b, 1]), (Subscript[c, 1]), (Subscript[d, 1]), (Subscript[p, 1]), (Subscript[s, 1])} with {(Subscript[a, 2]), (Subscript[b, 2]), (Subscript[c, 2]), (Subscript[d, 2]), (Subscript[p, 2]), (Subscript[s, 2])}.

(9) If (Subscript[s, 1]) = 0 goto step 2. If (Subscript[s, 1]) = 1, add intrv((Subscript[a, 1]), (Subscript[b, 1]), (Subscript[c, 1]), (Subscript[d, 1])) to rootlist, else put interval data {(Subscript[a, 1]), (Subscript[b, 1]), (Subscript[c, 1]), (Subscript[d, 1]), (Subscript[p, 1]), (Subscript[s, 1])} on intervalstack.

(10) If (Subscript[s, 2]) = 0 goto step 2. If (Subscript[s, 2]) = 1, add intrv((Subscript[a, 2]), (Subscript[b, 2]), (Subscript[c, 2]), (Subscript[d, 2])) to rootlist, else put interval data {(Subscript[a, 2]), (Subscript[b, 2]), (Subscript[c, 2]), (Subscript[d, 2]), (Subscript[p, 2]), (Subscript[s, 2])} on intervalstack. Go to step 2.

Correctness of the above algorithm was proved and its computational complexity was discussed in [Akr1978b, Akr1980, Akr1989].

Memory usage

Conjecture

sgc (p (x + 1)) + 
        sgc ((((x + 1))^n) p ((1/(x + 1)))) ≤ 
      sgc (p (x))

implies that the number of sign changes in any interval on intervalstack is at least equal to the total number of sign changes in all intervals above it. Since the number of sign changes in the top interval is at least 2, and the total number of sign changes in all intervals on stack is at most the degree n of f, the number of possible levels in
intervalstack is at most (Subscript[log, 2]) n. Therefore the maximal number of transformed polynomials we need to keep at any given time is at most 1 + (Subscript[log, 2]) n.

Experimental Results

We compare performance of our continued fraction (CF) algorithm, and the algorithm REL described in [RZ2002]. We have implemented both algorithms as a part of Mathematica kernel. They both use the same implementation of Shaw and Traub's algorithm for Taylor shifts (see [GG1997]). As benchmark examples we use Chebyshev, Laguerre, Wilkinson, and Mignotte polynomials used in [RZ2002], as well as three types of randomly generated polynomials used in [ABS1994].

All computations were done on a 850 MHz Athlon PC with 256 MB RAM. The memory used data was obtained using Mathematica MaxMemoryUsed command, so it includes the total memory used by Mathematica kernel. The startup size of Mathematica kernel is 1.6 MB.

The results given for random polynomials were averaged over sets of 5 random polynomials each, both methods were tested on the same sets of randomly generated polynomials.

Special polynomials

Here CF is faster by factors ranging from around 3 for Chebyshev polynomials to 50000 for Mignotte polynomials. The case of Mignotte polynomials is especially advantageous for our continued fractions method, because there is a point with a very simple continued fraction expansion (namely (1/5)), which lies between the two close roots. For Chebyshev polynomials we used the fact that the polynomials are even and so with both methods we isolated only the positive roots.

({{Polynomial, Degree,   Number  

          of Roots,          CF

          
          T (s)/
              M (MB), \(\(\ \ \ \ \ \ \ \ \ \)\(REL

          T \((s)\)/M \((MB)\)\)\)}, {Chebyshev, 1000, 1000, 2172/9.2, 7368/8.5}, {Chebyshev, 1200, 1200, 4851/12.8, 15660/11.8}, {Laguerre, 900, 900, 3790/8.7, 22169/14.1}, {Laguerre, 1000, 1000, 6210/10.4, 34024/17.1}, {Wilkinson, 800, 800, 73.4/3.24, 3244/10}, {Wilkinson, 900, 900, 143/3.66, 5402/12.5}, {Wilkinson, 1000, 1000, 256/4.1, 8284/15.1}, {Mignotte, 300, 4, 0.12/1.75, 803/7.7}, {Mignotte, 400, 4, 0.22/1.77, 3422/15.8}, {Mignotte, 600, 4, 0.54/1.89, 26245/49.1}})

Polynomials with randomly generated coefficients

Here CF was faster by factors between 1.5 and 5.

({{Coefficients

          (bit length), Degree,   Number  

          of Roots,          CF

          
          T (s)/
              M (MB), \(\(\ \ \ \ \ \ \ \ \ \)\(REL

          T \((s)\)/M \((MB)\)\)\)}, {10, 500, 3.6, 0.78/2.2, 1.66/2.81}, {10, 1000, 4.4, 6.67/3.75, 34.2/7.5}, {10, 2000, 5.6, 215/11.4, 562/22.8}, {1000, 500, 3.2, 0.56/2.28, 2.19/2.97}, {1000, 1000, 3.6, 12.7/5.1, 31.4/6.5}, {1000, 2000, 6.0, 329/14.2, 510/24.3}})

Monic polynomials with randomly generated coefficients

This case proved to be especially hard for REL. In this case CF was several thousand times faster. This is because such polynomials tend to have both very large and small roots, so an isolation method based on interval bisection starts with a very large interval, and needs to bisect it many times before it isolates the small roots. CF does not have this problem, because the size of its each ``step'' is based on an estimate of how far the next root is.

({{Coefficients

          (bit length), Degree,   Number  

          of Roots,          CF

          
          T (s)/
              M (MB), \(\(\ \ \ \ \ \ \ \ \ \)\(REL

          T \((s)\)/M \((MB)\)\)\)}, {10, 500, 5.2, 1.43/2.48, 8.84/3.84}, {10, 1000, 4.8, 7.12/3.74, 80.7/10.1}, {10, 2000, 6.8, 263/11.4, 1001/37.1}, {1000, 100, 4.4, 0.01/1.75, 56.8/5.5}, {1000, 200, 6.0, 0.086/1.93, 252/17}, {1000, 500, 5.6, 0.57/2.28, 1917/96.8}, {1000, 1000, 6.0, 25.5/5.2, >5000/?}})

Products of factors (x-randomly generated integer root)

Here CF was up to 25 times faster for small roots, but REL was up to 4 times faster for large roots. The latter being the only case when we found CF to be slower than REL.

({{       Roots

          (bit length), Degree,   Number  

          of Roots,          CF

          
          T (s)/
              M (MB), \(\(\ \ \ \ \ \ \ \ \ \)\(REL

          T \((s)\)/M \((MB)\)\)\)}, {10, 100, 100, 0.8/1.82, 0.61/1.92}, {10, 200, 200, 2.45/2.07, 10.1/2.64}, {10, 500, 500, 33.9/3.34, 878/8.4}, {1000, 20, 20, 0.12/1.88, 0.044/1.83}, {1000, 50, 50, 16.7/3.18, 4.27/2.86}, {1000, 100, 100, 550/8.9, 133/6.49}})

Conclusions

We have shown that our continued fraction root isolation algorithm is almost always faster than the algorithms based on interval bisection. Its bound on memory usage, given in terms of the number of transformed polynomials it needs to keep, is not much worse then for the algorithm REL presented in [RZ2002], and in practice its memory usage is often smaller than that of REL.

References

[Akr1978] Alkiviadis G. Akritas, A correction on a theorem by Uspensky, Bulletin of the Greek Mathematical Society,  19 (1978), pp. 278-285.
[Akr1978b] Alkiviadis G. Akritas, Vincent's theorem in algebraic manipulation, Ph.D. thesis, North Carolina State University, Raleigh, N.C., 1978.
[Akr1980] Alkiviadis G. Akritas, An implementation of Vincent's Theorem, Numerische Mathematik, \textbf{36}(1980), pp. 53-62.
[Akr1982] Alkiviadis G. Akritas, Reflections on a pair of theorems by Budan and Fourier, Mathematics Magazine, 55 (1982), pp. 292-298.
[Akr1984] Alkiviadis G. Akritas, Budan's theorem and its consequences, Sciences et Techniques en Perspective, 4 (1984), pp. 1-13.
[Akr1989] Alkiviadis G. Akritas, Elements of Computer Algebra-with Applications, Wiley, New York, NY, 1989. Available also in Russian, MIR Publishers, Moscow, 1994 (with new material).
[ABS1994] Alkiviadis G. Akritas, Alexei V. Bocharov and Adam W. Strzebonski, Implementation of real root isolation algorithms in Mathematica, Abstracts of the International Conference on Interval and Computer-Algebraic Methods in Science and Engineering (Interval'94), St.Petersburg, Russia, March 7-10, (1994), pp.23-27.
[AD1980] Alkiviadis G. Akritas and Stylianos D. Danielopoulos, On the complexity of algorithms for the translation of polynomials, Computing, 24 (1980), pp. 51-60.
[CA1976] George E. Collins and Alkiviadis G. Akritas, Polynomial real root isolation using Descartes' rule of signs, Proceedings of the 1976 ACM Symposium on Symbolic and Algebraic Computations, Yorktown Heights, N.Y., (1976), pp. 272-275.
[GG1997] Joachim von zur Gathen and Jurgen Gerhard, Fast Algorithms for Taylor Shifts and Certain Difference Equations, Proceedings of ISSAC'97, Maui, Hawaii, U.S.A., (1997), pp. 40-47.
[LT1972] S. Lang and H. Trotter, Continued fractions for some algebraic numbers, Journal fur die reine und angewandte Mathematik, 255 (1972), pp. 122-134.
[RZ2002] Fabrice Rouillier and Paul Zimmermann, Efficient isolation of polynomial real roots, preprint, 2002.
[U1948] J. V. Uspensky, Theory of Equations, McGraw-Hill, New York, 1948.
[V1836] A. J. H. Vincent, Sur la resolution des equations numeriques, Journal de Mathematiques Pures et Appliquees, 1 (1836), pp. 341-372.