RESEARCH ANNOUNCEMENT: FASTER FACTORIZATION OF TWO NUMBERS 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.

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.

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.

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.

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 at most n(lg n)^{1+o(1)}; the gcd takes time at most 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 at most 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}.

Date: 2004.10.09. 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.

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 at most m(lg m)^{2+o(1)} where m is the total length of g_{i−1}, g_{i}. The total time here is at most 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: Compute b mod g_{1}, b mod g_{2}, … with a scaled remainder tree; 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.

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.

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