# Brent Pollard ρ factorisation

% Remco Bloemen % 2014-02-27

In the previous post the Pollard ρ method was by far the fastest of the methods tested to factor an 64 bit integer. In this post I implement an improved version by John Pollard and Richard Brent. It is described in this paper:

Richard P. Brent (1980) – An Improved Monte Carlo Factorization Algorithm

The straightforward implementation in C++ is:

```
uint64 brent_pollard_factor(uint64 n)
{
const uint64 m = 1000;
uint64 a, x, y, ys, r, q, g;
do a = random() % n; while(a==0||a==n-2);
y = random() % n;
r = 1;
q = 1;
do
{
x = y;
for(uint64 i=0; i < r; i++)
{
// y = y² + a mod n
y = mul(y, y, n);
y += a;
if(y < a) y += (max_uint64 - n) + 1;
y %= n;
}
uint64 k =0;
do
{
for(uint64 i=0; i < m && i < r-k; i++)
{
ys = y;
// y = y² + a mod n
y = mul(y, y, n);
y += a;
if(y < a) y += (max_uint64 - n) + 1;
y %= n;
// q = q |x-y| mod n
q = mul(q, (x>y)?x-y:y-x, n);
}
g = gcd(q, n);
k += m;
} while(k < r && g == 1);
r <<= 1;
} while(g == 1);
if(g == n)
{
do
{
// ys = ys² + a mod n
ys = mul(ys, ys, n);
ys += a;
if(ys < a) ys += (max_uint64 - n) + 1;
ys %= n;
g = gcd((x>ys)?x-ys:ys-x, n);
} while(g == 1);
}
return g;
}
```

The parameter `m`

has been fixed to `1000`

, this seemed to be the
fastest choice.

Given a composite number `n`

it will return a non-trivial factor of `n`

.
To turn this method into a prime factorisation method this factor must
be further split up in factors until the numbers are prime. I develop
the function listed below to do this.

It maintains a stack of possibly composite factors and a list of prime
already discovered. It pops a possible composite from the stack, checks
its primality (with this method).
If the number is prime, it is removed from all the other factors on the
stack and put on the prime list, else, the number is split using either
trial division, if it is small, or `brent_pollard_factor`

.

```
vector<uint64> prime_factors(uint64 n)
{
vector<uint64> factors;
vector<uint64> primes;
uint64 factor = brent_pollard_factor(n);
factors.push_back(n / factor);
factors.push_back(factor);
do
{
uint64 m = factors[factors.size() - 1];
factors.pop_back();
if(m == 1) continue;
if(is_prime(m))
{
primes.push_back(m);
// Remove the prime from the other factors
for(int i=0; i < factors.size(); i++)
{
uint64 k = factors[i];
if(k % m == 0)
{
do k /= m; while(k % m == 0);
factors[i] = k;
}
}
}
else
{
factor = (m < 100) ? small_factor(m) : brent_pollard_factor(m);
factors.push_back(m / factor);
factors.push_back(factor);
}
} while(factors.size());
return primes;
}
```

## Benchmark and conclusion

Compared to the methods posted previously, this method can perform the
task in 0.3 s. A tenfold improvement over the previous best and three
thousand times faster than trial division. Also, the time spent in `gcd`

is less than 20%, much less than the previous 90%.