Another Pollard Rho Implementation

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;
}

Solution

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