We studied discrete logarithms in two previous exercises. Today we look at a third algorithm for computing discrete algorithms, invented by John Pollard in the mid 1970s. Our presentation follows that in the book *Prime Numbers: A Computational Perspective* by Richard Crandall and Carl Pomerance, which differs somewhat from other sources.

Our goal is to compute *l* (some browsers mess that up; it’s a lower-case ell, for “logarithm”) in the expression *g*^{l} ≡ *t* (mod *p*); here *p* is a prime greater than 3, *g* is an integer generator on the range 1 ≤ *g* < *p*, and *t* is an integer target on the range 1 ≤ *g* < *p*. Pollard takes a sequence of integer pairs (*a*_{i}, *b*_{i}) modulo (*p* − 1) and a sequence of integers *x*_{i} modulo *p* such that *x*_{i} = *t*^{ai} g^{bi} (mod *p*), beginning with *a*_{0} = *b*_{0} = 0 and *x*_{0} = 1. Then the rule for deriving the terms of the various sequences is:

- If 0 <
*x*_{i} < *p*/3, then *a*_{i+1} = (*a*_{i} + 1) mod (*p* − 1), *b*_{i+1} = *b*_{i}, and *x*_{i+1} = *t x*_{i} (mod *p*).
- If
*p*/3 < *x*_{i} < 2*p*/3, then *a*_{i+1} = 2 *a*_{i} mod (*p* − 1), *b*_{i+1} = 2 *b*_{i} mod (*p* − 1), and *x*_{i+1} = *x*_{i}^{2} mod *p*.
- If 2
*p*/3 < *x*_{i} < *p*, then *a*_{i+1} = *a*_{i}, *b*_{i+1} = (*b*_{i} + 1) mod (*p* − 1), and *x*_{i+1} = *g x*_{i} mod *p*.

Splitting the computation into three pieces “randomizes” the calculation, since the interval in which *x*_{i} is found has nothing to do with the logarithm. The sequences are computed until some *x*_{j} = *x*_{k}, at which point we have *t*^{aj} *g*^{bj} = *t*^{ak} *g*^{bk}. Then, if *a*_{j} − *a*_{j} is coprime to *p* − 1, we compute the discrete logarithm *l* as (*a*_{j} − *a*_{k}) *l* ≡ *b*_{k} − *b*_{j} (mod (*p* − 1)). However, if the greatest common divisor of *a*_{j} − *a*_{j} with *p* − 1 is *d* > 1, then we compute (*a*_{j} − *a*_{k}) *l*_{0} ≡ *b*_{k} − *b*_{j} (mod (*p* − 1) / *d*), and *l* = *l*_{0} + *m* (*p* − 1) / *d* for some *m* = 0, 1, …, *d* − 1, which must all be checked until the discrete logarithm is found.

Thus, Pollard’s rho algorithm consists of iterating the sequences until a match is found, for which we use Floyd’s cycle-finding algorithm, just as in Pollard’s rho algorithm for factoring integers. Here are outlines of the two algorithms, shown side-by-side to highlight the similarities:

# find d such that d | n # find l such that g**l = t (mod p)
function factor(n) function dlog(g, t, p)
func f(x) := (x*x+c) % n func f(x,a,b) := ... as above ...
t, h, d := 1, 1, 1 j := (1,0,0); k := f(1,0,0)
while d == 1 while j.x <> k.x
t = f(t) j(x,a,b) := f(j.x, j.a, j.b)
h = f(f(h)) k(x,a,b) := f(f(k.x, k.a, k.b))
d = gcd(t-h, n) d := gcd(j.a-k.a, p-1)
return d return l ... as above ...

Please pardon some abuse of notation; I hope the intent is clear. In the factoring algorithm, it is possible that *d* is the trivial factor *n*, in which case you must try again with a different constant in the *f* function; the logarithm function has no such possibility. Most of the time consumed in the computation is the modular multiplications in the calculations of the *x* sequence; the algorithm itself is O(sqrt *p*), the same as the baby-steps, giant-steps algorithm of a previous exercise, but the space requirement is only a small constant, rather than the O(sqrt *p*) space required of the previous algorithm. In practice, the random split is made into more than 3 pieces, which complicates the code but speeds the computation, as much as a 25% improvement on average.

Your task is to write a program that computes discrete logarithms using Pollard’s rho algorithm. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.