Pollard’s rho algorithm

In 1903 Frank Nelson Cole went to the blackboard at a meeting of a mathematical society, silently computed 2⁶⁷ − 1, then multiplied 193,707,721 by 761,838,257,287, and the results matched. He later admitted that finding the divisors had taken “three years of Sundays”. Trial division would have had to go up to 193,707,721, almost a hundred million divisions by odd numbers alone. Pollard’s rho algorithm finds the same divisor in fourteen thousand steps. It does not look for the divisor at all: it waits until a pseudo-random sequence gives it away by repeating modulo a prime unknown to us. Below it factors real numbers, from the Fermat number that Euler factored in 1732 to the 78-digit Fermat number F₈.

A step here is one evaluation of x ↦ x² + c modulo n, and gcd is how many times the greatest common divisor was computed. While the algorithm runs, the line under the answer shows how far it has “reached”: after s steps a divisor below (s/1.25)² would most likely have been found already. After a million steps that is 6·10¹¹, while trial division in a million divisions checks divisors only up to two million. On a 30-digit number with a 13-digit divisor rho spends about two million steps, on a 40-digit one with a 16-digit divisor about thirty million. The time of trial division in the same cases is hours and months. The browser factors the Fermat number F₈ with its 16-digit divisor in a few seconds, and with the trick described below several times faster.

The idea

Take the sequence x_i+1 = x_i² + c mod n and look at it modulo an unknown prime divisor p. There are only p possible values there, and by the birthday paradox a repeat x_i ≡ x_j (mod p) appears after about √p steps. It can be noticed without knowing p: at a repeat p divides x_i − x_j, and gcd(x_i − x_j, n) is greater than one.

A repeat modulo p comes much earlier than modulo n: there are n/p times as many values modulo n. So the gcd is usually p or a multiple of it, but not n itself. It remains not to try all pairs (i, j), whose number is the square of the number of steps, but to notice the repeat while keeping only a couple of numbers. For this there are cycle-finding techniques.

The letter ρ

Modulo a prime p the sequence is deterministic: the next value depends only on the current one. So after the first repeat it goes round in a circle. Draw its path, first a tail of μ values, then a cycle of length λ, and you get the Greek letter ρ; hence the name. Below are real sequences modulo 101, 1009 and 10,007. Change c and the starting value and watch the tortoise and the hare run.

Floyd. The tortoise takes one step, the hare two, and on every move their values are compared. While both are on the tail, the hare runs ahead, but once on the cycle it starts gaining on the tortoise by one point per move and meets it before completing a lap. In Pollard’s algorithm “compare” means computing gcd(x_i − x_2i, n): we cannot see equality modulo p, but we can see a common divisor.

Brent. The hare runs one step at a time, and the tortoise stands still and jumps to the hare when it has run 1, 2, 4, 8, … steps. Sooner or later the stretch between jumps becomes longer than the cycle, and the hare comes to the tortoise. Each value is computed once rather than three times, so there are fewer evaluations of f.

The histogram under the letter is the length of ρ for each of the p possible starts with the chosen c. The line is what a random mapping on p points predicts, with mean √(πp/2) ≈ 1.25√p. The mapping x² + c is not random at all, but behaves similarly, and the whole algorithm rests on this.

Why √p and not √n

The algorithm rests on three claims and one hypothesis. The first two claims are about the sequence itself, the third about truly random numbers, and the hypothesis is that x² + c behaves like a random one.

Hence the letter ρ: modulo p there are only p values, the next is determined by the current one, so within p steps a value repeats, and after that the sequence goes round in a circle. The tail μ and the cycle λ together are no longer than p.

It remains to estimate μ + λ, the index of the first repeat. If values were chosen randomly and independently from p possible ones, this would be the birthday problem.

At k = 1.2√p the right-hand side is below a half, at k = 3√p about one per cent: a repeat almost certainly happens within a few √p steps. The mean can be computed exactly; it is about √(πp/2) ≈ 1.25√p. And modulo n there are n/p times as many values, so a repeat there takes √(n/p) times longer, which is why the gcd almost always gives p and not n.

        p  tries  rho/sqrt p  Floyd/sqrt p  Brent/sqrt p
    10007   1000       1.224         3.042         1.927
   100003    600       1.299         3.236         2.049
  1000003    300       1.249         3.098         2.028
 10000019    100       1.318         3.187         2.105

Floyd needs about 3.1√p evaluations of f on average, one by the tortoise and two by the hare on each move, and Brent about 2.0√p, a third fewer. Brent himself in his 1980 paper estimated the gain at 36% in cycle finding and 24% in factorisation.

Code

Factorisation in full: small divisors and the primality of pieces are left to the Miller–Rabin test, and a composite piece is split by rho. The Brent variant has two improvements. The tortoise waits at powers of two, as in the scene. And the gcd is not computed at every step: the differences x − y are multiplied modulo n, and once in 128 steps one gcd of the product is taken. If p divides any of the differences, it divides the product too. The rare risk that all divisors fall into one batch at once and the gcd equals n is removed by repeating the last batch one step at a time.

2^64 + 1 = 274177 * 67280421310721
    Floyd: 2424 steps, 808 gcds; Brent: 1918 steps, 13 gcds; trial division: 4101232 divisions
2^67 - 1 = 193707721 * 761838257287
    Floyd: 16584 steps, 5528 gcds; Brent: 13822 steps, 50 gcds; trial division: 96853860 divisions
p * q = 999999937 * 1000000007
    Floyd: 36777 steps, 12259 gcds; Brent: 28670 steps, 102 gcds; trial division: 499999968 divisions

Brent takes about a fifth fewer steps, and over a hundred times fewer gcds: 50 instead of 5528 on 2⁶⁷ − 1. The gcd of two long numbers is a whole Euclidean algorithm, noticeably dearer than one multiplication, so the batches save more than the cycle finding itself. The number 2⁶⁴ + 1 = 274,177 · 67,280,421,310,721 was factored by Fortuné Landry in 1880, without any computing machines.

What it costs

The steps number about √p, where p is the smallest prime divisor, and each step takes a multiplication and a reduction modulo n: O(log² n) by the schoolbook method. In all O(√p · log² n), with O(1) memory. The worst case is two factors of the same length, p ≈ √n, and then the steps number about n^1/4. For comparison, trial division makes up to p divisions if p is to be found, and up to √n if primality is to be proved.

The measurement shows the main property. Numbers n of 64, 128 and 256 bits, each with a smallest divisor p of given length, and the time to the first divisor:

All three lines have a slope of about ½: a p sixteen million times larger costs several thousand times more. The length of n does not affect the slope; it only shifts the line up, because multiplying long numbers costs more: with p around 2³⁰ a 256-bit number takes 1.7 times longer than a 64-bit one in Python and almost 5 times longer in Go. Another surprise: Go with math/big is no faster than Python here. The numbers are short, and the time goes not on arithmetic but on calls and memory allocation; on uint64 with bits.Mul64 Go would beat Python several times over.

A trick for Fermat numbers

Sometimes something is known about the divisors in advance. Euler proved that the prime divisors of a Fermat number F_m = 2^2m + 1 have the form k·2^m+1 + 1, and Édouard Lucas in 1878 strengthened this to k·2^m+2 + 1. For a divisor of F₈ this means that p − 1 is divisible by 2¹⁰. Take the mapping x¹⁰²⁴ + c instead of x² + c. Modulo such a p the power x¹⁰²⁴ takes only about p/1024 distinct values, the sequence runs over a set a thousand times smaller, and the cycle is √1024 = 32 times shorter. A step costs ten squarings instead of one, but the gain is greater. In the first scene turn this trick off with the checkbox and compare: without it F₈ needs about thirty million steps, with it fewer than a million.

On exactly this number Richard Brent and John Pollard showed the power of the method in 1980. F₈ was known to be composite, but not a single divisor was known. A modification of the rho algorithm found 1,238,926,361,552,897 in two hours on a Univac 1100/42, about 23 million multiplications modulo F₈. The primality of the remaining 62-digit factor was first proved by Hugh Williams.

What makes factorisation hard

The rho algorithm is exactly as fast as the smallest divisor is small. A random large number almost always has small divisors, and rho together with the Miller–Rabin test factors it down to the last, usually prime, piece. Hard are numbers with no small divisors at all: a product of two primes of the same length. For a 30-digit product of two 15-digit primes rho needs about 6·10⁷ steps, seconds in the browser; for a 60-digit one about 10¹⁵, for a 600-digit one 10¹⁵⁰. These are exactly the numbers used in RSA, and more on strength is in the basic course chapter “A secret in plain sight” and in the section “Breaking in” of “Mathematics, the Queen of the Sciences”.

Other methods split the same task differently. The p − 1 method, also by Pollard (1974), finds p instantly if all prime divisors of p − 1 are small, and is powerless otherwise. Hendrik Lenstra’s elliptic curve method (1985) also “looks for a small divisor”, but its time grows more slowly than any power of p; it finds divisors of dozens of digits. The quadratic sieve and the number field sieve work differently: their time depends on the size of n, not p, and for numbers of hundreds of digits without small divisors they are unrivalled. Pollard was the first to propose the number field sieve too, in 1988; in 1990 it was used to factor the ninth Fermat number, of 155 digits. In practice the methods line up into a pipeline: trial division, rho, elliptic curves, sieve.

The same trick for logarithms

In 1978 Pollard applied the same idea to the discrete logarithm: find x with g^x = h in a group of N elements. A pseudo-random walk over the group again cycles in √N steps, and x is extracted from the coincidence of two points of the path. For elliptic curves no better general methods are known, which is why the strength of a curve over a field of 2²⁵⁶ elements is estimated at 2¹²⁸ steps; more in the chapter “Curves in your browser”.

Where it works

• The factor utility from GNU coreutils factors numbers just like the program above: trial division by a table of small primes, a primality check by Miller–Rabin and the Lucas test, then Pollard–Brent rho or the SQUFOF method for the remaining composite piece.
• PARI/GP: the factorint function combines the rho algorithm, Shanks’ SQUFOF, the elliptic curve method and the quadratic sieve.
• Programming contests: factoring numbers up to 10¹⁸ is a standard problem for the Miller–Rabin and rho combination with multiplication through a 128-bit product.
• Cryptanalysis of elliptic curves: rho for logarithms, parallelised across many machines, is the main way of record discrete logarithm computations on curves.

History

John Pollard published the method in 1975 in a short paper, “A Monte Carlo method for factorization”. The cycle in it was caught by a method that Donald Knuth in 1969 attributed to Robert Floyd without a reference; it does not appear in Floyd’s own work and may be folklore. In 1980 Richard Brent replaced it with his own cycle finding with the tortoise at powers of two, and the same year, together with Pollard, factored F₈. Divisors of Fermat and Mersenne numbers that used to cost “years of Sundays” became a matter of hours of machine time, and now of fractions of a second in the browser.

Problems

Restore the start of the Cunningham table: factor 2^k ± 1 for all k up to a hundred. Trial division up to a million is enough for most rows, but not for all.