An Improved Algorithm for Diophantine Equations in One Variable

Adam Strzebonski

We present a new algorithm for computing integer roots of univariate polynomials. For a polynomial f ∈ ℤ(t) we can find the set of its integer roots in a time polynomial in the size of the
sparse encoding off. An algorithm to do this was given by F. Cucker, P. Koiran, S. Smale in 1999. Their paper introduces a polynomial time method for computing sign of f at integral points, and it finds integer roots of f by isolating the real roots of all polynomials in the sparse derivative sequence of f, up to unit-length intervals. We propose to isolate the roots of f using a sparse variant  of Fourier`s theorem, and show that our algorithm requires a smaller number of polynomial sign evaluations. We present an empirical comparison of our implementations of the original and of the improved algorithm, and of an algorithm using modular root
finding and Hensel lifting, suitable for dense polynomials.

Introduction

Let

f = (Subscript[a, 1]) (Subscript[(t^e), 1]) + … + (Subscript[a, k]) 
          (Subscript[(t^e), k])

where (Subscript[e, 1]) > … > (Subscript[e, k]), and \(TraditionalForm\`a\_i ∈ ℤ\{0}\), and let us denote

minexp (f) := (Subscript[e, k])

The sparse derivative sequence of f is the sequence (Subscript[f, 1]), …, (Subscript[f, k])defined by

(Subscript[f, 1]) := f/(t^(minexp (f)))
(Subscript[f, i + 1]) := (Subscript[f, i])'/(t^(minexp ((Subscript[f, i])')))

The algorithm defined in [CKS] (F.Cucker, P.Koiran, S.Smale, ``A Polynomial Time Algorithm for Diophantine Equations in One Variable'', J.Symbolic Comp., 27 (1999), 21-29), finds integer roots of f by isolating the real roots of polynomials of the sparse derivative sequence of f, taken in the reverse order, up to intervals (u,  u + 1) or [u, u], with u ∈ ℤ. The method is based on the fact that, by Rolle's theorem, a polynomial can have at most one root between two subsequent roots of its derivative. In the worst case the algorithm may require isolating (k^2)  roots.

Our algorithm isolates the roots of f directly, using a ``sparse variant'' of Fourier's theorem proven in the next section. It requires isolating at most  2 k - 2 roots of f or polynomials in its sparse derivative sequence.

We will also present an algorithm finding roots of f modulo a prime and using Hensel lifting to obtain the integer roots. Since f may have (Subscript[e, 1]) roots modulo a prime p, even if  k = 2, for instance if f = (t^p) - t, this algorithm is exponential in the size of f. However, our experiments suggest that for dense polynomials this method is the fastest.

Root Isolation Theorem

Let (Subscript[f, 1]), …, (Subscript[f, k]) : 
        I⟶ℝ be a sequence of differentiable functions defined in an open interval I ⊂ ℝ, and such
that for all 1 ≤ i ≤ k - 1, and for all x ∈ I

sign ((Subscript[f, i + 1]) (x)) = 
      sign ((Subscript[f, i])' (x))

and

(Subscript[f, k]) (x) ≠ 0

For  x ∈ I, and 1 ≤ p ≤ q ≤ k let (Subscript[sgc, p, q])(x) denote the number of sign changes in the sequence (Subscript[f, p])(x), …, (Subscript[f, q])(x) with terms equal to zero removed, and let sgc(x) := (Subscript[sgc, 1, k])(x).

Theorem: For any a ∈ I

\(lim\+\(x -> \(a\^+\)\)\ sgc \((x)\) = 
      sgc \((a)\)\)
\(lim\+\(x -> \(a\^-\)\)\ sgc \((x)\) = 
      sgc \((a)\) + r + 2  s\)

where \(TraditionalForm\`r ∈ \(ℤ\_+\)\), \(TraditionalForm\`s ∈ \(ℤ\_+\) ⋃ {0}\), and

(Subscript[f, 1]) (a) = … = (Subscript[f, r]) (a) = 
          0
(Subscript[f, r + 1]) (a) ≠ 0

Moreover, s = 0, unless there is t > r + 1, such that (Subscript[f, t])(a) = 0.

Corollary: Let (Subscript[f, 1]), …, (Subscript[f, k]) be the sparse derivative sequence of a polynomial f, let 0 ≤ a < b, and let 1 ≤ p ≤ q ≤ k, be such that (Subscript[sgc, p, k])(a) = (Subscript[sgc, p, k])(b). Then

(1) (Subscript[f, p]) has a constant nonzero sign on (a, b]

(2)  (Subscript[sgc, p, k])(c) = (Subscript[sgc, p, k])(a) for any a ≤ c ≤ b

(3) (Subscript[sgc, 1, p])(a) - (Subscript[sgc, 1, p])(b) = 
        r + 2  s, where r is the number of roots of f in  (a, b], counted with multiplicities, and s is a nonnegative integer. In particular, if (Subscript[sgc, 1, p])(a) = (Subscript[sgc, 1, p])(b) then f has no roots in (a, b], and if  (Subscript[sgc, 1, p])(a) - (Subscript[sgc, 1, p])(b) = 1 then f has exactly one root in (a, b], and sign(f(a)) ≠ sign(f(b)).

Corollary: A univariate polynomial with k > 0 nonzero terms has at most k - 1positive roots and at most k - 1 negative roots, and hence at most 2 k - 1 real roots. For completeness, let us note that the polynomial

f = (∏_(i = 0)^(k - 1))((t^2) - 
            (i^2))

has exactly k nonzero terms and 2  k - 1 real roots.

The Algorithm

We present an algorithm computing the positive integer roots of a polynomial \(TraditionalForm\`f ∈ ℤ[t]\{0}\). To find the negative integer roots of f we compute the positive integer roots of f(-t). To compute the sign of a polynomial we use the algorithm described in [CKS]. As a positive integer root bound for f, we take the smaller of the bound used in [CKS]  (absolute value of the trailing coefficient of f), and the positive root bound described in J. R. Johnson, ``Algorithms for Polynomial Real Root Isolation'', in B. Caviness, J. Johnson (eds.), Quantifier Elimination and Cylindrical Algebraic Decomposition, Springer-Verlag 1998, 269-299. We also use the following two
algorithms.

BinarySearch(f,a,b): for \(TraditionalForm\`f ∈ ℤ[t]\{0}\) and a, b ∈ ℤ,  such that f has exactly one root in (a, b), and f(a) f(b) < 0, finds out whether the root is an integer, and if yes returns it. The algorithm uses binary search
with integer subdivision points, and stops when it finds an integer root, or isolates the root to an interval of length 1.

ExhaustiveSearch(f,a,b,r): for \(TraditionalForm\`f ∈ ℤ[t]\{0}\) and a, b ∈ ℤ, finds the integer roots of f  in [a, b], assuming that there is at most r of them. The algorithm computes the sign of f at subsequent integers in [a, b], and stops when it finds r roots, or when it has checked all integers in [a, b].

Algorithm 1: Given \(TraditionalForm\`f ∈ ℤ[t]\{0}\) returns the list of positive integer roots of f.

(1) Compute the sparse derivative sequence (Subscript[f, 1]), …, (Subscript[f, k]) of f, a positive integer root bound rb for (Subscript[f, 1]), and set rootlist to an empty list. While (Subscript[f, 1])(rb) = 0 and rb > 0 add rb to rootlist and decrement rb. Put interval (0, rb) on intervalstack.

(2) If intervalstack is empty, return rootlist, else take an interval (a, b) off intervalstack.

(3) If sgc(a) = sgc(b) there are no roots of f in (a, b). Go to step 2.

(4) If  sgc(a) - sgc(b) = 1, compute BinarySearch((Subscript[f, 1]), a, b). If it finds an integer root add it to rootlist. Go to step 2.

(5) If b - a ≤ k, compute ExhaustiveSearch((Subscript[f, 1]), a, b, sgc(a) - sgc(b)). Add the roots found to rootlist. Go to step 2.

(6) Choose integers a < c = d < b. While (Subscript[f, 1])(c) = 0 and c > a add c to rootlist and decrement c. While (Subscript[f, 1])(d) = 0
and d < b add d to rootlist and increment d. Put intervals (a, c) and (d, b) on intervalstack and go to step 2.

Remark: There are several efficiency improvents that can be made to the algorithm described above.

(a) sgc(0) is equal to the number of sign changes in the coefficient list of f. We can compute it before computing the sparse derivative sequence, and if it is 0 or 1, we know that there are no roots, or we can find the only root using BinarySearch. Similarly,if rb is small (rb < k in our implementation), we use ExhaustiveSearch on [1, rb], without computing the sparse derivative sequence.

(b) We can construct the loop somewhat differently, so that we check whether sgc(a) - sgc(b) is 0 or 1, or whether b - a is small before adding (a, b) to intervalstack. This way intervalstack contains only intervals for which  sgc(a) - sgc(b) > 1, and hence its height never exceeds (k - 1)/2.

(c) Keep the computed signs of (Subscript[f, 1]), …, (Subscript[f, k]) at interval's endpoints together with the interval, so that we don't need to recompute them.

(d) For every interval (a, b) which we need to subdivide, compute the smallest p satisfying the conditions of Corollary. Then for any point c subdividing (a, b) we only need to compute the signs of (Subscript[f, 1])(c), …, (Subscript[f, p - 1])(c).

(e) In point 6 of the algorithm we choose c = floor ((a + b)/2), except that since the changes in value of sgc tend to happen close to zero, for intervals with  a = 0we choose c to be closer to zero than b/2.

Complexity Analysis

Let us compute the number of polynomial sign evaluations required by our algorithm (let us call it FIR, for Fourier Integer Roots), and by the algorithm described in [CKS] (let us call it RIR for Rolle Integer Roots). Let \(TraditionalForm\`f ∈ ℤ[t]\{0}\), let k be the number of nonzero terms in f, and let  M be a bound for integer roots of polynomials in the sparse derivative sequence (Subscript[f, 1]), …, (Subscript[f, k]) of f. For simplicity we take a common root bound for all  (Subscript[f, 1]), …, (Subscript[f, k]). Since the worst case estimate for size of coefficients of polynomials (Subscript[f, i]) grows with i, taking a common estimate is ``to the advantage'' of RIR.

FIR isolates at most 2  k - 2 changes in value of sgc, from an interval of length of at most M to an interval of length
of at least 1, evaluating signs of at most k polynomials for each interval subdivision. Hence, the total number of polynomial sign evaluations is bounded by O((k^2) log(M)).

To isolate roots of (Subscript[f, k - m + 1]) RIR needs to evaluate (Subscript[f, k - m + 1]) at at most (((m - 1))^2) endpoints of intervals isolating roots
of (Subscript[f, k - m + 2]), …, (Subscript[f, k]), and then it needs to isolate at most 2  m - 1 roots of (Subscript[f, k - m + 1]) from an interval of length of at most M to intervals of length 1, evaluating 1 polynomial sign for each subdivision. Hence, the number of polynomial sign evaluations necessary for isolating roots of  (Subscript[f, k - m + 1]), assuming the roots of (Subscript[f, k - m + 2]), …, (Subscript[f, k])  have already been isolated, is bounded by O((m^2) + m log(M)). Therefore the total number of polynomial sign evaluations necessary for finding
the integer roots of f using RIR is bounded by O((k^3) + (k^2) log(M)).

Clearly when k ≫ log(M) FIR has a lower complexity than RIR. However, the complexity computed above is a worst case estimate. In practice we quite often get large intervals with sgc(a) - sgc(b) = 1, in which case in the further subdivision of the interval we need only to compute signs of  (Subscript[f, 1]). Also, the use of Remark helps to lower the number of polynomial signs that need to be computed with each interval subdivision.

The Dense Polynomial Case

Another possible approach to finding integer roots of a polynomial f is to find the roots of f modulo a prime p, and then use Hensel lifting to find the integer roots. Such an algorithm is exponential in the size of f because the number of roots modulo p may be as large as the degree of f. Therefore the algorithm is not suitable for sparse polynomials. However, it works quite well when f is dense.

Algorithm 2: Given \(TraditionalForm\`f ∈ ℤ[t]\{0}\) returns the list of  integer roots of f.

(1) Compute the square free part g of f, and an integer root bound rb for g.

(2) Select a prime p such that g is square free modulo p.

(3) Compute roots (Subscript[r, 1]), …, (Subscript[r, s]) of g modulo p, by computing h = gcd(g, (t^p) - t), and factoring it using the algorithm
described in [Z] R. Zippel, ``Effective Polynomial Computation'', Kluwer Academic Publishers, Boston/Dodrecht/London 1993, section 18.3.

(4) Lift (Subscript[r, 1]), …, (Subscript[r, s]) to roots (Subscript[r, 1])', …, (Subscript[r, s])' of g mod (p^((2^m))), with (p^((2^m))) ≥ 2  rb + 1, using the quadratic iteration Hensel lifting described in [Z], section 16.2.

(5) Return those of  (Subscript[r, 1])' - (p^((2^m))), …, 
      (Subscript[r, s]) - (p^((2^m)))', (Subscript[r, 1])', …, (Subscript[r, s])' which are in the interval [-rb, rb] and are roots of g.

The integer root bound is the greater of the positive and negative integer root bounds computed as in Algorithm 1. In our implementation whenever the root bound is less than 10 we find the integer roots using a direct search. In step 2 we try an increasing sequence of primes, until we find one for which g is square free modulo p. Primes for which g is not square free modulopdivide the discriminant of g therefore there is only a finite number of them.

Experimental Results

In this section we present an empirical comparison of our implementations of Algorithm 1 (FIR, for Fourier Integer Roots), Algorithm 2 (HIR, for Hensel Integer Roots), and of the algorithm described in [CKS] (RIR, for Rolle Integer Roots). Both FIR and RIR use the refinement proposed in section 4 of [CKS]. All three algorithms have been
implemented as a part of Mathematica kernel. The computations have been done on a 700 MHz Pentium III laptop computer with 256 MB of RAM. Time is given in seconds, columns marked st\# show the number of polynomial sign tests used in FIR and RIR. The results given are averages for sets of 10 randomly generated polynomials.

We have used the following three types of randomly generated polynomials.

Random sparse polynomials of degree n with k integer roots

\(SP \((n, k, 
          p)\) := \((\(a\_1\) \(t\^\(n - k\)\) \(\(+∑\)\+\(i = 
                    2\)\%\(p - 1\)\) \(a\_i\) t\^e\_i + 
            a\_p)\) \(∏\+\(j = 1\)\%k\((t - 
              r\_j)\)\)\)

where (Subscript[a, i]) and (Subscript[r, j]) are randomly generated 100-bit integers, and (Subscript[r, j]) are all distinct, and 0 < (Subscript[e, i]) < n - k are randomly generated distinct exponents.

For all random sparse polynomials we tried FIR was the fastest algorithm, except for polynomials with the highest density (29%and 43%), where HIR was faster. The speed advantage of FIR over RIR grows with the number of terms. For very high degree polynomials we see that the number of polynomial sign tests performed both by FIR and RIR decreases. This is a result of using the refinement proposed in section 4 of [CKS].

({{Polynomial, Deg, Terms, Roots, Time

          HIR, Time

          FIR, Time

          RIR, st  #

          FIR, st  #

          RIR}, {SP ((10^3), 1, 10), (10^3), 19.9, 1, 1.1, 0.005, 0.005, 110, 317}, {SP ((10^3), 1, 10), (10^3), 29.2, 2, 1.11, 0.011, 0.085, 274, 3616}, {SP ((10^3), 1, 10), (10^3), 49.2, 4, 1.14, 0.055, 0.96, 799, 27299}, {SP ((10^3), 1, 10), (10^3), 89.3, 8, 1.19, 0.22, 7.21, 2198, 122415}, {SP ((10^3), 1, 10), (10^3), 157.1, 16, 1.31, 0.88, 65.6, 5196, 575958}, {SP ((10^3), 1, 10), (10^3), 294.9, 32, 1.64, 4.72, >300, 12825, ?}, {SP ((10^2), 4, 10), (10^2), 43.1, 4, 0.028, 0.034, 0.68, 722, 22672}, {SP ((10^4), 4, 10), (10^4), 50, 4, 40.4, 0.33, 0.38, 456, 2874}, {SP ((10^5), 4, 10), (10^5), 50, 4, >300, 5.61, 5.65, 418, 1031}, {SP ((10^6), 4, 10), (10^6), 50, 4, >300, 76.5, 76.5, 420, 1031}})

Random dense polynomials of degree n with k integer roots

DP (n, 
          k) := ((∑_(i = 0)^(n - k))(Subscript[a, i]) 
              (t^i)) (∏_(j = 1)^k)(t - 
              (Subscript[r, j]))

where (Subscript[a, i]) and (Subscript[r, j]) are randomly generated 100-bit integers, and (Subscript[r, j]) are all distinct.

When the polynomials are dense, HIR is clearly the fastest method. FIR is slower, but still much faster than RIR.

({{Polynomial, Deg, Roots, Time

          HIR, Time

          FIR, Time

          RIR, st  #

          FIR, st  #

          RIR}, {DP (128, 128), 128, 128, 1.09, 28.3, >300, 22372, ?}, {DP (128, 32), 128, 32, 0.23, 1.33, 130, 6548, 734241}, {DP (128, 8), 128, 8, 0.086, 0.26, 12.2, 2337, 218165}, {DP (128, 2), 128, 2, 0.052, 0.13, 2.1, 879, 58697}, {DP (128, 0), 128, 0, 0.008, 0.015, 0.11, 82, 545}})

Sparse polynomials with k + 1 terms and the maximal number 2  k + 1of integer roots

PP (k, p) := 
      t (∏_(j = 1)^k)((t^p) - ((Subscript[r, j])^p))

where pis an even integer and (Subscript[r, j]) are randomly generated 10-bit integers.

For random sparse polynomials with the maximal number of roots the fastest algorithm is FIR. HIR gives the same timings as FIR for smaller degree polynomials, but for higher degrees it becomes much slower. RIR is slower than FIR by a factor which decreases with the degree.

({{Polynomial, Deg, Terms, Roots, Time

          HIR, Time

          FIR, Time

          RIR, st  #

          FIR, st  #

          RIR}, {PP (5, 2), (10^1), 6, 10, 0.003, 0.003, 0.007, 139, 337}, {PP (5, 20), (10^2), 6, 10, 0.006, 0.007, 0.012, 128, 314}, {PP (5, 200), (10^3), 6, 10, 0.018, 0.018, 0.034, 126, 268}, {PP (5, 2000), (10^4), 6, 10, 16.3, 0.46, 0.51, 129, 142}, {PP (5, 20000), (10^5), 6, 10, >300, 14.5, 16.3, 126, 142}})