Another Pollard Rho Implementation

Question

In an attempt to solve the 3rd problem on project Euler (https://projecteuler.net/problem=3), I decided to implement Pollard's Rho algorithm (at least part of it, I'm planning on including the cycling later). The odd thing is that it works for numbers such as: 82123(factor = 41) and 16843009(factor 257). However when I try the project Euler number: 600851475143, I end up getting 71 when the largest prime factor is 6857. Here's my implementation(sorry for wall of code and lack of type casting):

#include <iostream>
#include <math.h>
#include <vector>

using namespace std;

long long int gcd(long long int a,long long int b);
long long int f(long long int x);

int main()
{
long long int i, x, y, N, factor, iterations = 0, counter = 0;
vector<long long int>factors;

factor = 1;
x = 631;
N = 600851475143;
factors.push_back(x);

while (factor == 1)
{
    y = f(x);
    y = y % N;
    factors.push_back(y);
    cout << "\niteration" << iterations << ":\t";
    i = 0;
    while (factor == 1 && (i < factors.size() - 1))
    {
        factor = gcd(abs(factors.back() - factors[i]), N);
        cout << factor << " ";
        i++;
    }
    x = y;
    //factor = 2;
    iterations++;
}

system("PAUSE");
return 0;
}

long long int gcd(long long int a, long long int b)
{
    long long int remainder;
    do
    {
    remainder = a % b;
    a = b;
    b = remainder;
    } while (remainder != 0);
    return a;
}

long long int f(long long int x)
{
//x = x*x * 1024 + 32767;
x = x*x + 1;
return x;
}

Answer 1

Here is an implementation of the Pollard's Rho algorithm written in C :

#include <inttypes.h>
#include <stdio.h>

typedef unsigned long long int u64;

u64 xor_random(u64 *s) {
    // A shift-register generator has a reproducible behavior across platforms.
    return *s ^= *s << 13, *s ^= *s >> 7, *s ^= *s << 17;
}

u64 xor_rand(u64 *seed, u64 min, u64 max) {
    // Return a random number within the specified range.
    return min + xor_random(seed) % (max - min + 1);
}

u64 pollard_rho(const u64 n, u64 *seed) {
    // Factorize N using the given seed to explore different sequences.
    u64 gcd = 1, a, b, c, i = 0, j = 1, x = 1, y = xor_rand(seed, 1, n - 1);
    for (; gcd == 1; ++i) {
        if (i == j) {
            if (j >> 18)
                break; // Timeout.
            j <<= 1;
            x = y;
        }
        a = y, b = y;
        for (y = 0; a; a & 1 ? b >= n - y ? y -= n : 0, y += b : 0, a >>= 1, (c = b) >= n - b ? c -= n : 0, b += c);
        y = (1 + y) % n;
        for (a = y > x ? y - x : x - y, b = n; (a %= b) && (b %= a););
        gcd = a | b;
    }
    return gcd;
}

int main(void) {
    const u64 numbers[] = {2175399929, 3069110221, 3927376457, 4650101551, 5718902089, 7893908123, 12264592849, 11938853153, 16306295659, 18237736631, 27663517949, 24443522911, 56923600907, 40791521263, 58216633199, 90280393909, 112425186329, 93498091213, 165705820351, 181488359447, 190753785143, 485939193757, 546916314293, 600851475143, 783487523561, 1032054630209, 799528427851, 1634360587319, 2023847347003, 1211940773833, 2800852171889, 2379264132979, 4310259270017, 7062649623461,
                           4414944843157, 7217908164707, 10920663358949, 14861402852993, 13054694846939, 32587772045321, 34418268654697, 24789309740719, 43751511444127, 67534381283503, 41155200889573, 126067218469847, 138868141684771, 120891030787517, 280827709682249, 153620750352803, 257422857142049, 314272081997693, 358481446433269, 404834966358103, 1083562243979291, 1045283729555563, 942435566692003, 1485240510648113, 2079374807499893, 1790397676985989, 3699064157742493, 3088761540503339,
                           3880047121088473, 6712112064424927, 5616519031783393, 8424277609928153, 12711893637547799, 13518311715881923, 11739219502337287, 22884231005938619, 30299474156982179, 32099290931348149, 71001924011460361, 48678265289549869, 61134905784265699, 90613499469027191, 88722250082959777, 114164062836748159, 253250470462687693, 278113325659607021, 174486836002714597, 418568314297273043, 522903562044842057, 388028774959168009, 878780904690558347, 809039149885743529, 777488646266872559,
                           1286793546054645163, 2259679290083711291, 1942672432348730897, 4545497999217627281, 2872387085447613407, 3090892722614326129, 9130947667582191569, 7787096734992841589, 7183925161407098971, 16981999183835246593LLU, 16452702208736988367LLU, 12388934618705660233LLU};
    const int count = sizeof(numbers) / sizeof(*numbers);
    u64 seed = 1, factor;
    for (int i = 0; i < count; ++i) {
        do factor = pollard_rho(numbers[i], &seed);
        while (factor == 1 || factor == numbers[i]);
        printf("%20" PRIu64 " has a factor %-20" PRIu64 "\n", numbers[i], factor);
        if (factor == 1 || factor == numbers[i] || numbers[i] % factor != 0) {
            printf("Error\n");
            return 1;
        }
    }
}

You can call this algorithm recursively to end up with a complete factorization of N, but be aware that the Pollard's Rho does not stop when :

• given a prime number.
• given a perfect power.

This example is taken from a more complete solution on GitHub that includes :

• a deterministic Miller-Rabin primality checker.
• a perfect power checker.
• a trial division algorithm.

The source code is into the public domain.

Answer 2

Pollard's rho algorithm guarantees nothing. It doesn't guarantee to find the largest factor. It doesn't guarantee that any factor it finds is prime. It doesn't even guarantee to find a factor at all. The rho algorithm is probabilistic; it will probably find a factor, but not necessarily. Since your function returns a factor, it works.

That said, your implementation isn't very good. It's not necessary to store all previous values of the function, and compute the gcd to each every time through the loop. Here is pseudocode for a better version of the function:

function rho(n)
    for c from 1 to infinity
        h, t := 1, 1
        repeat
            h := (h*h+c) % n # the hare runs ...
            h := (h*h+c) % n # ... twice as fast
            t := (t*t+c) % n # as the tortoise
            g := gcd(t-h, n)
        while g == 1
        if g < n then return g

This function returns a single factor of n, which may be either prime or composite. It stores only two values of the random sequence, and stops when it finds a cycle (when g == n), restarting with a different random sequence (by incrementing c). Otherwise it keeps going until it finds a factor, which shouldn't take too long as long as you limit the input to 64-bit integers. Find more factors by applying rho to the remaining cofactor, or if the factor that is found is composite, stopping when all the prime factors have been found.

By the way, you don't need Pollard's rho algorithm to solve Project Euler #3; simple trial division is sufficient. This algorithm finds all the prime factors of a number, from which you can extract the largest:

function factors(n)
    f := 2
    while f * f <= n
        while n % f == 0
            print f
            n := n / f
        f := f + 1
    if n > 1 then print n