A note on Pollard’s ρ-method

Joan Boyar

1 Introduction

The modern factoring algorithms are all probabilistic. One of the first of these is due to John Pollard and is called Pollard’s ρ-method. See

John Pollard, A Monte Carlo method for factorization, BIT 15 (1975) 331-334.

Pollard’s ρ-method for factoring integers appears to have running time O(p^{1/2}), where p is the smallest factor of n, the number being factored. For composite numbers, this is O(n^{1/4}). The analysis of the method is heuristic and is based on the “birthday paradox”.

2 The birthday paradox

The birthday paradox expresses a fact which is surprising to most people when they first hear it.

If there are at least 23 randomly chosen people in a room, the probability that there are at least two with the same birthday is greater than 1/2.

In fact, as the number of people gets larger than 23, the probability of finding at least two people with the same birthday goes up quickly. For example, if there are 32 people, the probability is greater than 3/4.

Why does this happen? Let’s look at the probability P that no two people in the room have the same birthday. Assume that there are k people in the room, and assume that the probability that any particular person has a birthday on any given day is 1/365 (so we are assuming that all years have exactly 365 days). First, number the people arbitrarily, from 1 to k. If the people numbered 1 through i − 1 all have distinct birthdays, the probability that person i has a different birthday from any of these is \frac{365−(i−1)}{365}. Thus, the probability that all k people have different birthdays is P = ∏_{i=1}^{k} \frac{366−i}{365}.

More generally, we might have N possible values that could be assigned to k different objects. (In the case of birthdays, N = 365.) If we choose k objects at random and each object is equally likely to have any of these N possible values, what is the probability that none have the same value? The answer is Q(N, k) = ∏_{i=1}^{k} \frac{N+1−i}{N} =

∏_{j=0}^{k−1}(1 − \frac{j}{N}) ≤ ∏_{j=0}^{k−1} e^{−j/N} = e^{−k(k−1)/2N}. Thus, one can approximate Q(N, k) ≈ e^{− \frac{k^{2}}{2N}}. If we want to achieve some fixed probability q, the number of objects needed is approximately \sqrt{−2N · log q}. One can also show that the expected number of objects needed in order to get at least two with equal values is bounded by c\sqrt{N}, where, for large N, c tends to \sqrt{π/2} ≈ 1.2533.

3 The algorithm

Pollard’s ρ-method consists of producing integers modulo the n being factored until the algorithm finds two which are congruent modulo some factor p of n. Then computing z, the difference of these two, and computing the greatest common divisor of z and n will give p, or some multiple of p which is also a factor of n. If the numbers are produced are completely random, and if the first very first pair which are congruent modulo p is found, the birthday paradox gives us that one expects to have to produce O(\sqrt{p}) different integers modulo n. This is not exactly what happens in Pollard’s ρ-method, so the argument is only heuristic. The integers modulo n are produced using a seed x_{1}, chosen randomly from the numbers modulo n, and the pseudorandom number generator

x_{i+1} = x_{i}^{2} + 1 (mod n).

These numbers appear to be random enough so that the algorithm actually runs in the expected amount of time.

But how do we determine when we have found two values x_{i} and x_{j} which are congruent modulo p? To determine if two values x_{i} and x_{j} are congruent modulo p, we compute gcd(x_{j} − x_{i}, n) and check if the result is greater than one or not. If the result is greater than one and x_{i} ≠ x_{j}, then this will give a nontrivial factor of n. Of course, we may be so unfortunate that x_{i} = x_{j}, in which case we would start over with a new seed, but this is very unlikely, compared to the probability that they are simply congruent modulo a factor of n.

Thus, in order to find two values which are congruent modulo a factor p of n, one could conceivably compare each new value with all previous values, doing the calculation gcd(x_{j} − x_{i}, n) each time. However, we “expect” to calculate Ω(\sqrt{p}) values before finding the first pair congruent modulo p, so doing all of these comparisons involves Ω(p) computations, which is no better than checking the integers 2, 3, ... consecutively until one finds a divisor of n.

In order to do better, one can use Floyd’s cycle detection algorithm. Notice that if x_{i} and x_{j} are congruent modulo p, then x_{i+1} and x_{j+1} are also congruent modulo p, so x_{i+k} and x_{j+k} are congruent modulo p for all k ≥ 0. Thus, as the values x_{1}, x_{2}, x_{3}, ... are computed, the first pair which is congruent modulo p will occur at some point, and after that the values will be periodic. Suppose the first pair which is congruent is x_{s} and x_{t}. Those values before x_{s} are considered the tail and those between x_{s} and x_{t−1} inclusive are called the cycle. One can draw this to look like the Greek letter ρ; hence the name for the algorithm.

Computing this generator and considering the values modulo p is a special case of computing the values for a generator

f : {1, 2, ..., q} → {1, 2, ..., q}.

Any sequence produced using a generator of this form must eventually be periodic since there are only a finite number of values possible, so we must eventually find two which are equal. How do we find the period? There is a cycle detection algorithm due to Floyd which can be described as follows: Compare the i^{th} and 2i^{th} elements produced, incrementing i, until you find an i where they are equal. Suppose that the seed is x. This algorithm can be expressed as follows: y ← x i ← 1 z ← f(x) j ← 2 while y ≠ z do

y ← f(y)

z ← f(f(z))
i ← i + 1
j ← j + 2

end output (i, j) and (y, z)

Which i and j do we find? When is x_{k} first equal to x_{2k}? Suppose that the tail of our sequence has length s − 1 and the cycle has length t − s. Then for arbitrary indices a and b with b > a, x_{a} = x_{b} if t − s divides b − a and if a ≥ s. The first time that t − s divides 2m − m and m ≥ s is when m = (t − s)(1 + ⌊(s − 1)/(t − s)⌋) ≤ (t − s) + (s − 1) = t − 1. So, if this is applied to the generator in Pollard’s ρ-method, the number of comparison is O(\sqrt{p}), and the algorithm becomes:

result ← n while result = n do

choose a random x ∈ {0, 1, 2, ..., n − 1} y ← x z ← (x ∗ x) + 1 (mod n) while gcd(z − y, n) = 1 do

y ← (y ∗ y) + 1 (mod n)

u ← (z ∗ z) + 1 (mod n)
z ← (u ∗ u) + 1 (mod n)

end

result = gcd(z − y, n)

output result

If we try to apply Pollard’s ρ-method to factor the number n = 1081, we might start with x_{1} = 53. Then, we compute

x_{1} = 53 x_{2} = 648 gcd(x_{2} − x_{1}, n) = 1
x_{2} = 648 x_{4} = 520 gcd(x_{4} − x_{2}, n) = 1
x_{3} = 477 x_{6} = 101 gcd(x_{6} − x_{3}, n) = 47

In this case 1081/47 = 23, so the algorithm found the larger prime factor before the smaller. This is more unusual when there is a bigger difference between the sizes of the factors.