RESEARCH ANNOUNCEMENT: FASTER FACTORIZATION INTO COPRIMES

DANIEL J. BERNSTEIN

Abstract. This paper presents an algorithm that, given positive integers a, b, computes the natural coprime base for {a, b} in time n(lg n)^{2+o(1)}, where n is the number of input bits. This paper also presents an algorithm that, given a set S of positive integers, computes the natural coprime base for S in time n(lg n)^{4+o(1)}.

1. Introduction

My previous paper [1] introduced an algorithm that, given a set S of positive integers, computes the natural coprime base cb S in time n(lg n)^{O(1)}, where n is the number of input bits. I made no attempt in [1] to optimize the exponent of lg n.

Section 2 of this paper presents an algorithm that computes cb{a, b} in time n(lg n)^{2+o(1)}. It is reasonable to conjecture that the limiting exponent 2 is optimal (for, e.g., a multitape Turing machine): one has cb{a, b} = {a, b} − {1} if and only if a, b are coprime; the well-known problem of checking coprimality has been stuck at n(lg n)^{2+o(1)} for thirty years.

Section 4 of this paper presents an algorithm that computes cb({a}∪Q), where Q is any coprime set, in time n(lg n)^{2+o(1)}. Substitute this algorithm into Algorithms 17.3 and 18.1 of [1] to obtain cb(P ∪ Q) in time n(lg n)^{3+o(1)} and cb S in time n(lg n)^{4+o(1)}. I’m not willing to conjecture that the 3 and 4 are optimal.

This is a very early draft. I’m confident in the basic structure of the algorithms, but there could be some silly omissions, and of course the proofs need vastly more detail.

2. Computing a coprime base for two positive integers

The following algorithm computes cb{a, b}, given positive integers a and b, in time n(lg n)^{2+o(1)}.

Step 1. Swap a, b if necessary so that a ≥ b. The algorithm will later reduce the input length by at least one third of the length of a. If a = 1, stop.

Date: 2004.11.03. Permanent ID of this document: 53a2e278e21bcbb7287b81c563995925. This is a preliminary version meant to announce ideas; it will be replaced by a final version meant to record the ideas for posterity. There may be big changes before the final version. Future readers should not be forced to look at preliminary versions, unless they want to check historical credits; if you cite a preliminary version, please repeat all ideas that you are using from it, so that the reader can skip it.

2000 Mathematics Subject Classification. Primary 11Y16.

The author was supported by the National Science Foundation under grant DMS–0140542, and by the Alfred P. Sloan Foundation.

Step 2. Compute a_{0} = a, g_{0} = gcd{a_{0}, b}, a_{1} = a_{0}/g_{0}, g_{1} = gcd{a_{1}, g_{0}^{2}}, a_{2} = a_{1}/g_{1}, g_{2} = gcd{a_{2}, g_{1}^{2}}, and so on, until g_{k} = 1.

For example, if a = 2^{100}3^{100} and b = 2^{137}3^{13}, compute a_{0} = 2^{100}3^{100}, g_{0} = 2^{100}3^{13}, a_{1} = 3^{87}, g_{1} = 3^{26}, a_{2} = 3^{61}, g_{2} = 3^{52}, a_{3} = 3^{9}, g_{3} = 3^{9}, a_{4} = 1, g_{4} = 1.

Lower level: The gcd inputs a_{i}, g_{i−1}^{2} are often highly unbalanced. To compute gcd{a_{i}, g_{i−1}^{2}}, first divide a_{i} by g_{i−1}^{2}, and then use any standard fast gcd algorithm to compute gcd{g_{i−1}^{2}, a_{i} mod g_{i−1}^{2}}. The division takes time n(lg n)^{1+o(1)}; the gcd takes time m(lg m)^{2+o(1)} where m is the length of g_{i−1}^{2}.

All the g’s together have length O(n), and k is at most about lg n, so the total time here is n(lg n)^{2+o(1)}.

Step 3. Compute x_{0} = g_{0}/gcd{g_{0}, g_{1}^{∞}}, x_{1} = g_{1}/gcd{g_{1}, g_{2}^{∞}}, and so on.

For example, if a = 2^{100}3^{100} and b = 2^{137}3^{13}, compute x_{0} = 2^{100}, x_{1} = 1, x_{2} = 1, x_{3} = 3^{9}.

Lower level: Compute each gcd{g_{i−1}, g_{i}^{∞}} as gcd{g_{i−1}, g_{i}^{2^{e_{i}}} mod g_{i−1}} where e_{i} is the smallest nonnegative integer satisfying 2^{2^{e_{i}}} ≥ g_{i−1}. The repeated squarings and gcd take time m(lg m)^{2+o(1)} where m is the total length of g_{i−1}, g_{i}. The total time here is n(lg n)^{2+o(1)}.

Step 4. Compute y_{0} = gcd{b, x_{0}^{∞}}, y_{1} = gcd{g_{0}, x_{1}^{∞}}, y_{2} = gcd{b, g_{1}, x_{2}^{∞}}, y_{3} = gcd{b, g_{2}, x_{3}^{∞}}, and so on.

For example, if a = 2^{100}3^{100} and b = 2^{137}3^{13}, the algorithm computes y_{0} = 2^{137}, y_{1} = 1, y_{2} = 1, y_{3} = 3^{13}.

Lower level: Use a scaled remainder tree to compute b mod g_{1}, b mod g_{2}, …; this takes time n(lg n)^{2+o(1)} since b, g_{1}, g_{2}, … together have length O(n). Then compute gcd{b, g_{1}} as gcd{b mod g_{1}, g_{1}}; compute gcd{b, g_{2}} as gcd{b mod g_{2}, g_{2}}; and so on.

Step 5. Recursively print cb{x_{0}, y_{0}/x_{0}}; cb{x_{1}, y_{1}}; cb{x_{2}, y_{2}}; and so on. Also print cb{a′} = {a′} − {1} and cb{b′} = {b′} − {1} where a′ = a/gcd{a, b^{∞}} and b′ = b/gcd{b, a^{∞}}. Note that a′ has already been computed; it equals a_{k}.

For example, if a = 2^{100}3^{100} and b = 2^{137}3^{13}, recursively print cb{2^{100}, 2^{37}} = {2} and cb{3^{9}, 3^{13}} = {3}. Also print cb{1} = {} and cb{1} = {}. The complete output is {2, 3}.

I claim that x_{0}y_{0}, x_{1}y_{1}, …, a′, b′ are coprime; that a = a′x_{0}x_{1}y_{1}x_{2}y_{2}^{3}x_{3}y_{3}^{7} ⋯; that b = b′y_{0}y_{1}y_{2}y_{3} ⋯; and that y_{0}x_{1}y_{1}x_{2}y_{2} ⋯, the product of inputs to the recursive calls, is at most ab/a^{1/3} ≤ (ab)^{5/6}. Each of these facts can be checked from the following table of ord_{p} values, expressed in terms of e = ord_{p} a and f = ord_{p} b:

g_{0} g_{1} g_{2}  g_{3} … x_{0} y_{0} x_{1} y_{1} x_{2} y_{2} … a′ b′
0 0 0  0 … 0 0 0 0 0  0 … 0 f  if e = 0
e 0 0  0 … e f 0 0 0  0 … 0 0  if 0 < e ≤ f
f e−f 0  0 … 0 0 e−f f 0  0 … 0 0  if f < e ≤ 3f
f 2f e−3f 0 … 0 0 0 0 e−3f f … 0 0  if 3f < e ≤ 7f
⋮
0 0 0  0 … 0 0 0 0 0  0 … e 0  if f = 0 < e

Consequently the outputs of the algorithm are coprime; a and b are products of powers of the outputs; and the recursion multiplies the total time by a bounded factor.

Note that one can easily factor a, b over cb{a, b} by tracing the factorizations a = a′x_{0}x_{1}y_{1}x_{2}y_{2}^{3}x_{3}y_{3}^{7} ⋯ and b = b′y_{0}y_{1}y_{2}y_{3} ⋯ through the recursion.

3. An algorithm without a catchy name

The following algorithm computes gcd{s, p^{∞}} for each s in a multiset S and for each p in a nonempty coprime set P. It takes time (k + 1)n(lg n)^{2+o(1)} if #P ≤ 2^{k}.

See [2] and [3] for similar algorithms.

Step 1. If #P = 1: Find p ∈ P. Use a scaled remainder tree to compute p mod s for each s ∈ S. Compute gcd{s, p^{∞}} as gcd{s, (p mod s)^{∞}}. This takes time n(lg n)^{2+o(1)}.

Assume from now on that #P ≥ 2.

Step 2. Select Q ⊆ P with #Q = ⌊#P/2⌋. Use a product tree to compute y = ∏_{p∈Q} p. This takes time n(lg n)^{2+o(1)}.

Step 3. Use a scaled remainder tree to compute y mod s for each s ∈ S. This takes time n(lg n)^{2+o(1)}.

Step 4. Compute gcd{s, (y mod s)^{∞}} = gcd{s, y^{∞}} for each s ∈ S. This takes time n(lg n)^{2+o(1)}.

Step 5. Apply the algorithm recursively to {gcd{s, y^{∞}} : s ∈ S} and Q; separately handle each s for which gcd{s, y^{∞}} = 1. This produces gcd{gcd{s, y^{∞}}, p^{∞}} = gcd{s, p^{∞}} for each p ∈ Q.

Apply the algorithm recursively to {s/gcd{s, y^{∞}} : s ∈ S} and P −Q; separately handle each s for which s/gcd{s, y^{∞}} = 1. This produces gcd{s/gcd{s, y^{∞}}, p^{∞}} = gcd{s, p^{∞}} for each p ∈ P − Q.

The product of inputs at each level of recursion is exactly the original product of inputs, so the total input size at each level of recursion is O(n).

4. Extending a coprime base

The following algorithm computes cb({a} ∪ Q), where Q is coprime, in time n(lg n)^{2+o(1)}.

Step 1. Use a product tree to compute b = ∏_{q∈Q} q. This takes time n(lg n)^{2+o(1)}.

Step 2. Define, and compute, a_{0}, g_{0}, a_{1}, g_{1}, …, x_{0}, x_{1}, …, y_{0}, y_{1}, …, a′ exactly as in Section 2. This takes time n(lg n)^{2+o(1)}.

Step 3. Write y_{i}(q) for gcd{q, y_{i}^{∞}}. Compute y_{0}(q), y_{1}(q), y_{2}(q), … for each q ∈ Q as explained in Section 3. This takes time n(lg n)^{2+o(1)} lg lg n = n(lg n)^{2+o(1)}. Check for and discard 1’s so that they do not slow down subsequent computations.

Step 4. Use product trees to compute z_{0} = ∏_{q} y_{0}(q), z_{1} = ∏_{q} y_{1}(q), etc. This takes time n(lg n)^{2+o(1)}.

Notice that z_{1}z_{2}^{3}z_{3}^{7} ⋯ divides a. Indeed, take any prime p dividing z_{1}z_{2}^{3}z_{3}^{7} · · ·. Recall that y_{1}, y_{2}, … are coprime, so p divides z_{i} for a unique i. Write e = ord_{p} a and f = ord_{p} b; then (2^{i} − 1)f < e ≤ (2^{i+1} − 1)f. Furthermore, p divides a unique q ∈ Q, and f = ord_{p} q = ord_{p} z_{i} by definition of b and z_{i}, so (2^{i}−1) ord_{p} z_{i} < ord_{p} a.

Step 5. Use a scaled remainder tree to compute a mod z_{0}, a mod z_{1}^{3}, a mod z_{2}^{7}, . . .. This takes time n(lg n)^{2+o(1)}, since the product z_{1}^{3}z_{2}^{7} ⋯ divides a^{3}.

Step 6. Use scaled remainder trees to compute (a mod z_{0}) mod y_{0}(q) = a mod y_{0}(q) for each q; (a mod z_{1}^{3}) mod y_{1}(q)^{3} = a mod y_{1}(q)^{3} for each q; (a mod z_{2}^{7}) mod y_{2}(q)^{7} = a mod y_{2}(q)^{7} for each q; etc. This takes time n(lg n)^{2+o(1)}.

Step 7. Compute gcd{a, y_{0}(q)}, gcd{a, y_{1}(q)^{3}}, gcd{a, y_{2}(q)^{7}}, etc. This takes time n(lg n)^{2+o(1)}.

Observe that gcd{a, y_{0}(q)}, gcd{a, y_{1}(q)^{3}}, gcd{a, y_{2}(q)^{7}}, etc. are the same as gcd{a, y_{0}(q)^{∞}}, gcd{a, y_{1}(q)^{∞}}, gcd{a, y_{2}(q)^{∞}}, etc. Consider, for example, a prime p dividing gcd{a, y_{2}(q)^{∞}}. Recall that 3f < e ≤ 7f where e = ord_{p} a and f = ord_{p} b = ord_{p} y_{2}(q); thus ord_{p} gcd{a, y_{2}(q)^{∞}} = e = ord_{p} gcd{a, y_{2}(q)^{7}}.

Step 8. Print cb{y_{i}(q), gcd{a, y_{i}(q)^{∞}}} for each i and q. Also print q/y_{0}(q)y_{1}(q) ⋯ for each q, and print a′. This takes time n(lg n)^{2+o(1)}.

References

[1] Daniel J. Bernstein, Factoring into coprimes in essentially linear time, to appear, Journal of

   Algorithms. ISSN 0196–6774. URL: http://cr.yp.to/papers.html. ID f32943f0bb67a9317d4
   021513f9eee5a.

[2] Daniel J. Bernstein, How to find small factors of integers, accepted to Mathematics of Computation; now being revamped. URL: http://cr.yp.to/papers.html.

[3] Daniel J. Bernstein, How to find smooth parts of integers, draft. URL: http://cr.yp.to/

   papers.html#smoothparts. ID 201a045d5bb24f43f0bd0d97fcf5355a.

Department of Mathematics, Statistics, and Computer Science (M/C 249), The University of Illinois at Chicago, Chicago, IL 60607–7045

E-mail address: djb@cr.yp.to