Let r be a positive integer. In this paper we study the problem of finding all r-power divisors of a given positive integer N , i.e., all positive integers p such that p r | N . Throughout the paper we write lg x := log 2 x, and unless otherwise specified, the "running time" of an algorithm refers to the number of bit operations it performs, or more formally, the number of steps executed by a deterministic multitape Turing machine [12]. We always assume the use of fast (quasilinear time) algorithms for basic integer arithmetic, i.e., for multiplication, division and GCD (see for example [19] or [4]).
Our main result is the following theorem.
Theorem 1.1 There is an explicit deterministic algorithm with the following properties. It takes as input an integer N 2 and a positive integer r lg N. Its output is a list of all positive integers p such that p r | N . Its running time is
Note that whenever we write in a complexity bound, we mean that the bound holds for all > 0, where the implied big-O constant may depend on .
The integers p referred to in Theorem 1.1 need not be prime. Of course, if p is a composite integer found by the algorithm, then the algorithm will incidentally determine the complete factorisation of p, as the prime divisors of p must also satisfy r | N .
The hypothesis r lg N does not really limit the applicability of the theorem: if r > lg N then the problem is trivial, as the only possible r-power divisor is 1.
The case r = 1 corresponds to the ordinary factoring problem, and in this case our algorithm is essentially equivalent to Coppersmith's method. As mentioned above, the complexity is O(N 1/4+ ), which does not improve on known results; currently, the fastest known deterministic factoring method has complexity O(N 1/5+ ) [9]. (It is interesting to ask whether the ideas behind [9] can be used to improve Theorem 1.1 when r 2. Our inquiries in this direction have been so far unsuccessful.)
In fact, when r = 1, Theorem 1.1 gives the more precise complexity bound O(N 1/4 (lg N ) 10+ ). It is apparently well known that Coppersmith's method has complexity O(N 1/4 (lg N ) C ) for some constant C > 0, but to the best of our knowledge, this is the first time in the literature that a particular value of C has been specified. On the other hand, we have not tried particularly hard to optimise the value of C, and it is likely that it can be improved. (One possible improvement is outlined in Remark 3.6.)
The authors of [1] were mainly interested in cryptographic applications, and this led them to focus on the case that N = p r q where p and q are roughly the same size. In this setting, they show that their algorithm is faster than ECM when r ≈ (log p) 1/2 , and that it even runs in polynomial time when r is as large as log p.
In this paper we take a different point of view: our goal is to determine the worst-case complexity, without any assumptions on the size of p, q or r.
To illustrate what difference this makes, consider again the case r = 2. This case is mentioned briefly in Section 6 of [1]. The authors point out that if N = p 2 q, where p and q are known to be about the same size, i.e., both p and q are within a constant factor of N 1/3 , then the running time of their method is O(N 1/9+ ), i.e., the number of search intervals is O(N 1/9+ ). However, in our more general setup, this is not the worst case. Rather, the worst case running time is O(N 1/8+ ), which occurs when searching for p ∼ N 1/4 and q ∼ N 1/2 .
More generally, for r 1 the worst case running time of O(N 1/4r+ ) stated in Theorem 1.1 occurs when p ∼ N 1/2r and q ∼ N 1/2 . By contrast, in the "balanced" situation considered in [1], where p, q ∼ N 1/(r+1) , one can show that the running time is only O(N 1/(r+1) 2 + ) (see Remark 3.5, and take θ = r/(r + 1)).
Although the core of our algorithm is essentially the same as the BDHG algorithm, our more general perspective requires us to make a few changes to their presentation. For instance, we cannot take the lattice dimension to be d ≈ r 2 (as is done in the main theorem of [1]), because this choice is suboptimal when r is small and fixed. Additional analysis is required to deal with potentially small values of p and q, and in general we must take more care than [1] in estimating certain quantities throughout the argument. For these reasons, we decided to give a self-contained presentation, not relying on the results in [1]. Remark 1.2 After this paper was accepted for publication, Dan Bernstein mentioned to us (personal communication) that the proof of Theorem 5.2 of [2] likely includes much of the argument needed to obtain our Theorem 1.1.
An important component of our algorithm, and of all algorithms pursuing Coppersmith's strategy, is a subroutine for finding all integer roots of a polynomial with integer coeffi-cients. This problem has received extensive attention in the literature, but we were unable to locate a clear statement of a deterministic complexity bound suitable for our purposes. For completeness, in Appendix A we give a detailed proof of the following result. For a polynomial f ∈ Z[x], we write f ∞ for the maximum of the absolute values of the coefficients of f . Theorem 1.3 Let b n 1 be integers. Given as input a polynomial f ∈ Z[x] of degree n such that f ∞ 2 b , we may find all of the integer roots of f in time
Note that this complexity bound is much stronger than what is needed for the application in this paper. However, it is still not quasilinear in the size of the input, which is O(nb). For further discussion, see Remarks A.8 and A.10.
In this section we recall the strategy of [1] for finding all integers p in a prescribed interval P -H p P + H such that p r | N , provided that H is not too large. We will prove the following theorem.
we may find a nonzero vector w in the lattice L := span Z (v 0 , . . . , v d-1 ) such that w 2 (d-1)/4 (det L)
Then in time
we may find a nonzero polynomial
Proof Set fi (y) := f i (Hy) ∈ Z[y] d , and let v i ∈ Z d be the vector whose j-th entry (for j = 0, . . . , d -1) is the coefficient of y j in fi (y). We will apply Lemma 2.2 to the vectors v 0 , . . . , v d-1 . Let
We claim that v i B for all i. First consider the case 0 i < rm. For any j = 0, . . . , i, the coefficient of y j in fi (y) = N m-i/r (P + Hy) i is equal to
, where we have used the hypotheses (2.3) and (2.2). For the case rm i < d, a similar argument shows that every coefficient of fi (y) = (P + Hy) i is bounded above by 2 d N d/r . Therefore every v i has coordinates bounded by 2 d N d/r+1 , and we conclude that v i B for all i.
Next we calculate the determinant of the lattice L := span Z (v 0 , . . . , v d-1 ), or equivalently, the determinant of the d × d integer matrix whose rows are given by the v i . Since deg fi (y) = i, this is a lower triangular matrix whose diagonal entries are given by the leading coefficients of the fi (y), namely
The determinant is the product of these leading coefficients, i.e.,
Invoking Lemma 2.2, we may compute a nonzero vector w ∈ L such that
in time O(d 5+ (lg B) 2+ ). Note that this time bound certainly dominates the cost of computing the vectors v i themselves, as the fi (y) may be computed by starting with f0 (y) = N m , and then successively multiplying by P + Hy and occasionally dividing by N . The hypotheses (2.1) and (2.2) imply that
The vector w corresponds to a nonzero polynomial h(y
in the Z-span of the fi (y). Applying the Cauchy-Schwartz inequality to the vectors w = ( h0 , . . . , hd-1 ) and (1, . . . , 1) yields
Moreover, each hj is divisible by H j , so we obtain in turn a polynomial h(x)
hj /H j for each j, the estimate (2.6) follows immediately.
Next we show that any r-power divisor that is sufficiently close to P corresponds to a root of h(x).
Proof We claim that h(x 0 ) is divisible by p rm . Since h(x) is a Z-linear combination of the f i (x) (where f i (x) is defined as in Proposition 2.3), it is enough to prove that p rm | f i (x 0 ) for all i. For the case 0 i < rm, we have f i (x 0 ) = N m-i/r p i . Since p r | N , we have p r(m-i/r ) p i | f i (x 0 ), and this implies that p rm | f i (x 0 ) because r i/r i. For the case i rm we have simply f i (x 0 ) = p i , which is certainly divisible by p rm .
On the other hand, the assumption -H x 0 H together with (2.6) and (2.4) implies that
We may now complete the proof of the main theorem of this section.
We first invoke Proposition 2.3, with inputs N , r, m, d, P and H, to find a polynomial h(x) satisfying (2.6). According to Proposition 2.4, we may then construct a list of candidates for p by finding all integer roots of h(x), which we do via Theorem 1.3.
To estimate the complexity of the root-finding step, recall from the proof of Proposition 2.4 that rm , so certainly |h j | < (P -H) rm for all j, and we obtain
Therefore in Theorem 1.3 we may take n := d and b := lg(N d/r ) = d r lg N . Note that the hypothesis b n is satisfied due to (2.1). The root-finding complexity is thus
which is negligible compared to the main bound (2.5). Finally, we must check each candidate for p to ensure that p r | N , which again requires negligible time.
We now consider the problem of searching for all integers p such that p r | N in an interval, say T p T , that is too large to be handled by a single application of Theorem 2.1. Given N , r, T and T , our strategy will be to choose parameters d, m and H, and then apply Theorem 2.1 to a sequence of subintervals of the form P -H p P + H that cover the target interval T p T . The overall running time will depend mainly on the number of subintervals, so our goal is to make H as large as possible. On the other hand, to ensure that the hypothesis (2.4) of Theorem 2.1 is satisfied, we also require that H < H where
The key issue is therefore to choose d and m to maximise H. For large d and m, the magnitude of H depends more or less on the ratio m/d; in fact, one finds that H is maximised when m/d ≈ lg T / lg N . The following result gives a simple formula for m (as a function of d) that is close to the optimal choice, and a corresponding explicit lower bound for H. and let H be defined as in (3.1). Then 1) , where
It is easy to check that d 1/(d-1) 2 1/2 < 3 for all d 2, so we find that 1) .
Continuing to estimate the exponent in this inequality, we obtain
where the last line follows from the inequalities 0 δ < 1 d-1 and 0 θ 1.
We may now estimate the time required to search a given interval T p T .
There is an explicit deterministic algorithm with the following properties. It takes as input positive integers N , r, T and T such that (2.1) holds (i.e., r lg N ) and such that 4
Its output is a list of all integers p in the interval T p T such that p r | N . Its running time is
where θ is defined as in (3.3). ). Also, the assumption T N 1/r implies that m (d -1)/r d/r, so (2.2) holds.
Let H be defined as in (3.1). Since
Moreover, (3.4) implies that θ 2 r/ lg N , so we have N θ 2 /r N 4/ lg N = 16 and hence H > 16/6 > 2.
Let H be the largest integer less than H, i.e., H := H -1. Then 2 H < H, and moreover, since H > 2, we also have
We may compute H by first approximating the d(d -1)-th root of the rational number
so this can all be done in time O((lg N ) 3+ ).
We now apply Theorem 2.1 with the parameters N , r, d, m, H, and with successively P = T + H, P = T + 3H, and so on, stopping when the interval [T, T ] has been exhausted by the subintervals [P -H, P + H]. The hypotheses (2.1), (2.2) and (2.3) have already been checked above, and (2.4) follows from our choice of H < H because P -H T . The number of subintervals is at most
Finally, since d lg N , the cost of each invocation of Theorem 2.1 is
Remark 3.3 A slightly better choice for d is to take d ≈ θ lg N , but this complicates the analysis and only improves the main result by a constant factor.
Finally we may prove the main theorem. Recall that we are given as input positive integers N 2 and r lg N , and we wish to find all positive integers p such that p r | N . Such divisors p must clearly lie in [1, N 1/r ].
We first check all p = 2, 3, . . . , 2 k by brute force, i.e., testing directly whether p r | N . Note that k may certainly be computed in time O((lg N ) 1+ ). To estimate the cost of checking up to 2 k , observe that k 2 (lg N )/r + 1 + 1.
Let C > 0 be an absolute constant such that 2 √ x + 1 + 1 x/4 + C for all x 1; it follows that k (lg N )/4r + C, and hence that 2 k N 1/4r . The cost of checking up to 2 k is therefore O(N 1/4r (lg N ) 1+ ), which is negligible compared to (1.1).
We now apply Proposition 3.2 to the intervals [2 k , 2 k+1 ], [2 k+1 , 2 k+2 ], and so on until we reach N 1/r , taking the last interval to be [2 j , N 1/r ] for suitable j. Since k 2 (lg N )/r, the precondition (3.4) is satisfied. For each interval we have (T -T )/T = O(1), and since θ ∈ [0, 1] we have
Therefore the cost of searching each interval is
Finally, the number of intervals is at most lg(N 1/r ) = O( 1 r lg N ).
For f, g ∈ Z[x], let res(f, g) ∈ Z denote the resultant of f and g.
x] be nonzero polynomials, and let n := deg f , m := deg g. Then
Proof (where the sum is taken over primes).
Proof Let ϑ(X) := p X log p denote the usual Chebyshev weighted prime counting function. The claim is that ϑ(X)/X > 1 3 log 2 (≈ 0.231) for all X 2. For X 101 this follows from [15, Thm. 10], which states that ϑ(X)/X > 0.84 for X 101. For 2 X < 101 the claim may be checked directly, for example by inspecting the graph of ϑ(X)/X in the reader's favourite computer algebra system.
The core idea of Loos' algorithm is the following well-known p-adic Hensel lifting strategy.
Proposition A. 4 Let p be a prime, let k be a positive integer, and let f ∈ (Z/p k Z)[x] be a polynomial of degree n 1. Let u ∈ {0, . . . , p -1}, and suppose that f (u) ≡ 0 (mod p), f (u) ≡ 0 (mod p).
Then there exists a unique v ∈ {0, . . . , p k -1} such that v ≡ u (mod p), f(v) ≡ 0 (mod p k ).
Proof We argue by induction on k. If k = 1, we simply take v = u. Now assume that k 2 and set := k/2 < k. By induction there exists a unique w ∈ {0, . . . , p -1} such that w ≡ u (mod p) and f (w) ≡ 0 (mod p ). We first establish uniqueness of v. Suppose that v has the desired properties, i.e., v ≡ u (mod p) and f (v) ≡ 0 (mod p k ). By the uniqueness of w, we must have v ≡ w (mod p ), say v = w + p t for some t ∈ {0, . . . , p k--1}. Expanding f around w, we find that
for some g ∈ (Z/p k Z)[x]. Substituting x = p t, and using the fact that p 2 ≡ 0 (mod p k ), we deduce that 0 ≡ f (w) + p tf (w) (mod p k ). Since f (w) ≡ f (u) ≡ 0 (mod p), we may solve for t to obtain
This establishes uniqueness of t (mod p k-), and hence of v (mod p k ). Moreover, the same calculation gives an explicit formula for v, proving existence.
To prove the complexity bound, suppose that we have already computed w and that we wish to lift to v. We first apply Horner's rule to compute f (w) and f (w) using O(n) arithmetic operations in Z/p k Z. Each such operation requires time O(lg(p k ) 1+ ). Similarly, we may invert f (w), and hence compute t and v, in time O(lg(p k ) 1+ ). Therefore, the time required to deduce v from w is O(n lg(p k ) 1+ ). The contributions from subsequent recursion levels form a geometric series, so the total cost of computing v from u is also O(n lg(p k ) 1+ ).
The next result shows how to find a reasonably small prime p for which the p-adic lifting strategy is guaranteed to succeed. Proof Since f is squarefree, the resultant D := res(f, f ) is nonzero. Our goal is to find a prime p such that p D.
First, by Lemma A.1 we have
One easily checks that (n + 1) n-1 n n for all n 1. Since f ∞ n f ∞ , we obtain
On the other hand, let X := 6bn + 6n lg n. If D is divisible by all primes p X, then D is divisible by their product, so lg |D| p X lg p > X/3 = 2nb + 2n lg n by Lemma A.3. This contradicts (A.1), so we conclude that there must exist a prime p X such that p D. To actually find such a prime, we run the following algorithm.
Step 1 (list primes). Make a list of all primes p Y for Y := 6nb + 6n lg n = O(nb). Using the sieve of Eratosthenes, this requires time
Step 2 (reduce f modulo primes). Let f p ∈ (Z/pZ)[x] denote the reduction of f modulo p. We compute f p for all p Y by applying a fast simultaneous modular reduction algorithm [19,Thm. 10.24] (i.e., using a remainder tree) to each coefficient of f . The bit size of the product of the primes is
, and the number of primes is certainly O(nb), so the cost per coefficient is O(n 1+ b 1+ ). The total cost over all coefficients is therefore O(n 2+ b 1+ ).
Step 3 (compute GCDs).
For each p Y , we compute gcd(f p , f p ) ∈ (Z/pZ)[x] using a quasilinear time GCD algorithm [19,Cor. 11 Finally, we return the least prime p for which f p = 0 and gcd(f p , f p ) = 1. As shown above, such a prime exists and satisfies p X.
We now give a deterministic root-finding algorithm for the squarefree case.
Proposition A.6 Let b n 1 be integers, and let f ∈ Z[x] be a squarefree polynomial of degree n such that f ∞ 2 b . Then we may find all integer roots of f in time
Proof As above, let f p ∈ (Z/pZ)[x] denote the reduction of f modulo p. We first invoke Proposition A.5 to find a prime p = O(nb) such that f p is nonzero and squarefree. Then we perform the following steps.
Step 1 (find roots mod p).
Compute the roots of f p in Z/pZ by brute force, i.e., by evaluating f p (i) for i = 0, . . . , p -1. Note that the integer roots of f correspond to distinct roots of f p , thanks to the squarefreeness of f p . Each f p (i) may be evaluated in time O(n 1+ b ), so the cost of this step is O(pn
Step 2 (find roots mod p k ).
Let f ∈ (Z/p k Z)[x] be the reduction of f modulo p k , where k is chosen to be the smallest integer such that
Applying Proposition A.4 to f , we lift each of the roots of f p found in Step 1 to a root of f . The uniqueness claim in Proposition A.4 implies that the resulting set of lifted roots in Z/p k Z includes the reductions modulo p k of all of the actual integer roots of f . To estimate the complexity, observe that p k p(n + 1) 1/2 2 n+b+1 , so
The cost of lifting each root is therefore O(n lg(p k ) 1+ ) = O(nb 1+ ), and the total cost of this step is O(n 2 b 1+ ).
For each root r ∈ Z/p k Z of f found in Step 2, we determine whether it arises from a genuine integer root of f as follows. We first lift r to a candidate root r * ∈ Z satisfying r * ≡ r (mod p k ) and r * ∈ [- 1 2 p k , 1 2 p k ). We next divide f by xr to obtain a polynomial ḡ ∈ (Z/p k Z)[x] such that f (x) = (xr)ḡ(x), and we lift ḡ to a polynomial g * ∈ Z[x] satisfying g * ≡ ḡ (mod p k ) and whose coefficients also all lie in [- 1 2 p k , 1 2 p k ). We then multiply xr * by g * (x) (in Z[x]) and check whether we obtain f . If so, then f (r * ) = 0, so r * must be the integer root corresponding to r. Otherwise, as we will see in the next paragraph, this r does not correspond to any integer root and we may ignore it. This procedure requires O(n) operations on integers of O(b) bits, i.e., O(nb 1+ ) bit operations, so the total cost over all roots is O(n 2 b 1+ ).
We now prove that the procedure described above does in fact find all integer roots. (The following argument is adapted from [19, §15.6].) Let r ∈ Z be a root of f . Then |r| f ∞ 2 b (as r divides the constant term of f ), and f factors as f (x) = (xr)g(x) for some g ∈ Z[x] satisfying g ∞ (n + 1) 1/2 2 n+b (by Lemma A.2). In particular, (A.2) ensures that |r| < p k /2 and g ∞ < p k /2. Let r ∈ Z/p k Z be the root of f corresponding to r, and let r * ∈ Z and g * ∈ Z[x] be the quantities computed in Step 3 for this r. Then r * ≡ r ≡ r (mod p k ), so we must have r * = r, since both sides lie in [- 1 2 p k , 1 2 p k ). Similarly, we have
so again we must have g * = g as the coefficients on both sides lie in [- 1 2 p k , 1 2 p k ). Therefore (xr * )g * (x) = f (x), and the procedure does indeed recover r.
Remark A. 7 Loos [11] imposes the additional requirement that p should not divide the leading coefficient of f , to ensure that deg f p = deg f . This is because he is searching for rational roots, not just integral roots. Our algorithm may also be easily adapted to this case.
Remark A.8 An interesting question is whether the complexity bound in Proposition A.6 can be improved to quasilinear, i.e., to O(n 1+ b 1+ ) bit operations. There are two main obstructions to this.
First, although f p ∈ (Z/pZ)[x] is squarefree for almost all primes p, it is difficult to predict in advance for which p this will occur. Consequently, in the proof of Proposition A.5 we were forced to test every prime up to O(nb). If we allow probabilistic algorithms, then we can find a suitable prime with high probability by randomly selecting p in the range 2 p X for some X = O(nb). This allows us to find a suitable p in expected quasilinear time. The complexity of Steps 1 and 2 in Proposition A.6, i.e., finding the roots modulo p and lifting them to Z/p k Z, can also be improved to (deterministic) quasilinear time by means of fast multipoint evaluation techniques. The resulting algorithm is quite similar to the root finding algorithm presented in [19,Thm. 15.21].
The second obstruction concerns Step 3 of Proposition A.6, namely, checking which of the candidate integer roots are in fact roots of f . We do not know how to carry out this step rigorously in quasilinear time, even allowing randomised algorithms. A similar issue occurs in [19,Thm. 15.21], where the last term of the given complexity bound corresponds in our notation to O(n 2+ b 1+ ). In the discussion following that theorem, von zur Gathen and Gerhard suggest testing the candidate roots modulo a small prime (different to p) as a way to quickly rule out incorrect candidates. This idea can be turned into a "Monte Carlo" algorithm: one would randomly choose a small prime q, compute f (r * ) (mod q) for all candidate roots r * , and declare the ones for which f (r * ) ≡ 0 (mod q) to be the true roots. We suspect that in this way one can obtain a quasilinear expected running time with an exponentially small probability of failure, but we have not checked the details.
In order to prove Theorem 1.3, we must first discuss the computation of GCDs in Z[x].
Let f, g ∈ Z[x] and let h := gcd(f, g). The idea of the "heuristic GCD" algorithm [5] is to use an integer GCD algorithm to compute gcd(f (N ), g(N )) for some choice of evaluation point N ∈ Z. If we are lucky, then gcd(f (N ), g(N )) will actually be equal to h(N ), and we may simply read off the coefficients of h(x) from h(N ), provided that N is not too small. However, it is possible for gcd(f (N ), g(N )) to contain extraneous factors unrelated to h(x). Usually these extraneous factors are small but in rare circumstances they can be very large. The algorithm can be made to tolerate extraneous factors up to a given size by taking larger values of N , at the expense of running more slowly. In practice, one usually takes a fairly small value of N , accepting a small chance of failure in order to get a fast algorithm. In the next result we work at the other extreme, taking N so large that the algorithm is guaranteed to work in all cases. We thereby obtain a GCD algorithm that is deterministic and completely rigorous (although unfortunately quite slow in practice).
. Assume that at least one of f and g is primitive. Define
and given f and g as input, we may compute h, f and g in time
Proof The inequalities (A.3) follow immediately from Lemma A.2, as f and g are divisors of f and g respectively, and h is a divisor of both.
For any integer N we have f
where δ(N ) := gcd( f (N ), g(N )).
We may bound the quantities δ(N )h i independently of N as follows. (This argument is adapted from [7,Thm. 4].) Since f and g are relatively prime, their resultant R := res( f , g) is nonzero, and there exist polynomials r, s ∈ Z
Therefore the quantities δ(N )h i are bounded by
We may now describe the actual algorithm for computing h, f and g.
Step 1. Set N := 2 c where c := (2n + 1) lg(n + 1) + 2n 2 + 2nb + n + b + 2.
As shown above, |δ(N )h i | 2 c-2 for all i. Notice also that c = O(nb). Step 2. We compute f (N ) and g(N ), which amounts to concatenating the coefficients of f and g with appropriate zero-padding (or one-padding in the case of negative coefficients). The integers f (N ) and g(N ) have bit size O(nc) = O(n 2 b), and the concatenation may be performed in linear time, i.e., in time O(n 2 b). Step 3. We compute gcd(f (N ), g(N )) using a quasilinear time GCD algorithm (see for example [18]). This requires time O((n 2 b) 1+ ) = O(n 2+ b 1+ ).
Step 4. We read off the coefficients δ(N )h i from (A.4). This is possible thanks to the bound |δ(N )h i | 2 c-2 , i.e., the coefficients do not "overlap". (In more detail, we may first read off δ(N )h 0 from the lowest c bits, i.e., by reading (A.4) modulo N = 2 c . After subtracting off this term, we may read off δ(N )h 1 from the next c bits, and so on.) This requires linear time O(n 2 b).
Step 5. Since we assumed that at least one of f and g is primitive, h is also primitive. We may therefore compute δ(N ) by taking the GCDs of the integers δ(N )h 0 , . . . , δ(N )h n . Each pairwise GCD requires time O(c 1+ ), so the total time required for this step is O(nc 1+ ) = O(n 2+ b 1+ ). Step 6. We now recover h(N ) = gcd(f (N ), g(N ))/δ(N ), and then f (N ) = f (N )/h(N ) and g(N ) = g(N )/h(N ). Using a quasilinear time integer division algorithm, this requires time O(n 2+ b 1+ ). Finally, we read off the coefficients of f and g from f (N ) and g(N ), in a similar manner to Step 4.
Remark A.10 To the best of the authors' knowledge, it is not known how to improve the complexity bound in Proposition A.9 to quasilinear without giving up on determinism. Several randomised quasilinear-time algorithms are known. Schönhage [16] analyses a variant of the heuristic GCD algorithm in which the evaluation point is chosen randomly. Another approach is to compute the GCD modulo a collection of randomly chosen small primes [19,Alg. 6.38].
We may now prove our main root-finding result.