Pollard's rho algorithm implementation

Question

I used std::multiset to determine how many times that factor appears. I'd like to know if it can be improved; I've focused on the algorithm rather than the code.

https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm

std::multiset<std::int64_t> pollards_rho(std::int64_t n)
{
 std::multiset<std::int64_t> factors{};
 if (n == 0 || n == 1)
 {
 return factors;
 }
 if (n < 0)
 {
 factors.insert(-1);
 n = -n;
 }
 // If n is even, add 2 and get rid of all even factors.
 if (n % 2 == 0)
 {
 factors.insert(2);
 while (n % 2 == 0)
 {
 n /= 2;
 }
 if (n == 1)
 {
 return factors;
 }
 }
 // Select an x_0 uniformly at random from [2, n - 1].
 // x_0 is chosen as 2 for simplicity.
 //
 // Floyd's cycle-finding algorithm.
 // x => x_i
 // y => x_i+1
 std::int64_t x{2}, y{x}, d{1};
 // f(x_i+1) = x_i^2 + 1 mod n
 auto f = [n](const std::int64_t num) noexcept
 {
 return (num * num + 1) % n;
 };
 // If we obtain a factor greater than 1, we are done.
 while (d == 1)
 {
 x = f(x);
 y = f(f(y));
 d = std::gcd(std::abs(x - y), n);
 }
 /**
 * if (d == n)
 * {
 * return std::nullopt;
 * }
 * else
 * {
 * return d;
 * }
 */
 // The algorithm normally returns either a non-trivial factor d or failure,
 // but we need to find all prime factors.
 if (d == n)
 {
 // If d and n are the same, we have a prime number.
 factors.insert(n);
 }
 else
 {
 // Look for other factors.
 factors.merge(pollards_rho(d));
 factors.merge(pollards_rho(n / d));
 }
 return factors;
}

Answer 1

I used std::multiset to determine how many times that factor appears

It's an option. You could also use a plain old map to store a factor (key) and its exponent (value):

• In case we only care about factors and not their exponents, the exponents are easier to ignore that way.
• In case we do care about the exponents of factors, we just get them, instead of having to count how many times a factor appears in a multiset.
• "High" exponents don't cost as much duplication, but really an exponent can never be high, at worst you could have 63 copies of 2 in your multiset.
• a map is a less rare/unusual data structure than multiset.

// If n is even, add 2 and get rid of all even factors.
if (n % 2 == 0)
{
 factors.insert(2);
 while (n % 2 == 0)
 {
 n /= 2;
 }
 if (n == 1)
 {
 return factors;
 }
}

This does not determine how many times the factor of 2 appears, it's only added once. By the way with C++20 you could use std::countr_zero instead of a loop, that combines well with returning a map instead of a multiset.

// f(x_i+1) = x_i^2 + 1 mod n
auto f = [n](const std::int64_t num) noexcept
{
 return (num * num + 1) % n;
};

This does not handle the case when num * num overflows, which it can easily do. You could fix this by restricting the input to fit in an int32_t (but still do this arithmetic in f with int64_t - the square of a 32-bit integer fits in 64 bits so it works), or implement a proper 64-bit modular multiplication, that's a bit tricky but possible.

// The algorithm normally returns either a non-trivial factor d or failure,
// but we need to find all prime factors.

Careful: Pollard's Rho can also spuriously fail to find a factor when its input is not prime.

// If d and n are the same, we have a prime number.

So this is a bug. You can get d == n when n is actually composite, but your code would decline to factor it.

• ozan: "You can get d == n when n is actually composite, but your code would decline to factor it." Is this the limitation of the algorithm or did I implemented it badly?
• user555045: @ozan the algorithm inherently works that way
• ozan: I see, I also fixed adding factors of 2 you mentioned, looking for the other issues. Thanks.