C++ - Brent-Pollard rho algorithm infinite loop

Question

I have this following function. I got it from two sites and tried to adapt it into my own but it didn't work very well.

When I test unsigned long int max - 2 or to but it as a number 4294967293 it puts the following code into an infinite loop where ys will keep coming back to the same value.

I am very new to factorisation algorithms and step by step help into why I am getting the infinite loop would be great.

The following code is just my "rho" function. I have another function called gdc which is the same as every other gdc recursive function out there.

unsigned long int rho(unsigned long int n) 
{
    if (n == 1) return n;
    if (n % 2 == 0) return 2;

    unsigned long int y = rand() % n;
    unsigned long int x;
    unsigned long long int ys = y;
    unsigned long int c;
    do
        c = rand() % n;
    while (c == 0 || c == n - 2);
    unsigned long int m = 1000;
    unsigned long int d = 1;
    unsigned long int q = 1;
    unsigned long int r = 1;
    while (d == 1)
    {
        x = y;
        for (int i = 0; i < r; i++)
        {
            y = y * y % n;
            y += c;
            if (y < c)
                y += (std::numeric_limits<unsigned long>::max() - n) + 1;
            y %= n;
        }
        int j = 0;
        while (j < r && d == 1)
        {
            ys = y;
            for (int i = 0; i < m && i < (r-j); i++)
            {
                y = y * y % n;
                y += c;
                if (y < c)
                    y += (std::numeric_limits<unsigned long>::max() - n) + 1;
                y %= n;
                q *= ((x>y) ? x - y : y - x) % n;
            }
            d = gcd(q, n);
            j += m;
        }
        r *= 2;
    }
    if (d == n)
    {
        do
        {
            ys = ys * ys % n;
            std::cout << ys << std::endl;
            ys += c;
            if (ys < c)
                ys += (std::numeric_limits<unsigned long>::max() - n) + 1;
            ys %= n;
            d = gcd( ((x>ys) ? x - ys : ys - x) , n);
        } while (d == 1);
    }
    return d;
}

The examples I adapted my code from :

• Brent-Pollard ρ factorisation
• Pollard Rho Brent Integer Factorization

EDIT

I have done as Amd has suggested and tidied up my code and moved repetitive lines into helper functions. However I am still getting an infinite loop for the d==n part near the bottom. For some reason, f(ys) ends up essentially returning the same thing that it had returned earlier so it keeps cycling through a series of values.

uint64_t rho(uint64_t n)
{
    if (n == 1) return n;
    if (n % 2 == 0) return 2;

    uint64_t y = rand() % n;
    uint64_t x;
    unsigned long long int ys = y;
    uint64_t c;
    do c = rand() % n; while (c == 0 || c == n - 2);
    uint64_t m = 1000;
    uint64_t d = 1;
    uint64_t q = 1;
    uint64_t r = 1;
    do 
    {
        x = y;
        for (int i = 0; i <= r; i++)
            y = f(y, c, n);

        int j = 0;
        do 
        {
            ys = y;
            for (int i = 0; i <= min(m, r-j); i++)
            {
                y = f(y, c, n);
                q *= (abs(x,y) % n);
            }
            d = gcd(q, n);
            j += m;
        } while (j < r && d == 1);
        r *= 2;
    } while (d == 1);
    if (d == n)
    {
        do
        {
            ys = f(ys, c, n);
            d = gcd(abs(x, ys), n);
        } while (d == 1);
    }
    return d;
}

Answer 1

It should always terminate. By the time you get to that point in the algorithm, x will be within the cycle (since y started at x and went all the way around to detect the cycle using Brent's cycle detection). The value ys started at y which started at x, so it will also continue around the cycle and eventually meet up with x again (see Floyd or Tortoise and Hare cycle detection). At this point you will have gcd(ys,x)==x, and that final loop will terminate.

There are several bugs in the posted implementation which I believe may be causing the issue:

x = y;
for (int i = 0; i < r; i++)                // should be strictly less than
    ...

    ys = y;
    for (int i = 0; i < min(m, r-j); i++)  // again, strictly less than
    {
        y = f(y, c, n);
        q = (q*abs(x,y)) % n;  // needs "mod" operator AFTER multiplication
    }
    ...

You can also replace the initialization of c with

uint64_t c = (rand() % (n-3))+1

if you want a number in the range [1,n-3].

Here's the original paper by Richard P. Brent: An Improved Monte-Carlo Factorization Algorithm

Answer 2

Edit:
this works for me, and not trapped inside infinite loop:

#include<iostream>
#include<stdint.h>

#define min(a,b) (a<b?a:b)
#define abs(x,y) (x > y? x - y : y - x)

uint64_t gcd(uint64_t m, uint64_t n) {
    while (true) {
        int r = m % n;
        if (r == 0) {
            return n;
        }
        m = n;
        n = r;
    }
}
uint64_t f(uint64_t y, uint64_t c, uint64_t n) {
    y = (y * y) % n;
    y += c;
    if (y < c)
        y += (std::numeric_limits<uint32_t>::max() - n) + 1;
    y %= n;
    return y;
}

uint64_t rho(uint64_t n)
{
    if (n == 1) return n;
    if (n % 2 == 0) return 2;

    uint64_t y = rand() % n;
    uint64_t x;
    uint64_t ys = y;
    uint64_t c;
    do c = rand() % n; while (c == 0 || c == n - 2);
    uint64_t m = 1000;
    uint64_t d = 1;
    uint64_t q = 1;
    uint64_t r = 1;
    do
    {
        x = y;
        for (int i = 0; i <= r; i++)
            y = f(y, c, n);

        int j = 0;
        do
        {
            ys = y;
            for (int i = 0; i <= min(m, r - j); i++)
            {
                y = f(y, c, n);
                q *= (abs(x, y) % n);
            }
            d = gcd(q, n);
            j += m;
        } while (j < r && d == 1);
        r *= 2;
    } while (d == 1);
    if (d == n)
    {
        do
        {
            ys = f(ys, c, n);
            d = gcd(abs(x, ys), n);
        } while (d == 1);
    }
    return d;
}
int main() {
    std::cout << rho(std::numeric_limits<uint32_t>::max() - 2) << "\n";     //9241
}

Old:
according to an improved monte carlo factorization algorithm: http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf
pp:182-183
error: ys = x * x % n;
correct: ys = ys * ys % n; // ys=f(ys)

please use uint32_t or uint64_t ... from stdint.h for good and clean coding style:

#include<iostream>
#include<stdint.h>

int gcd(int m, int n) {
    while (true) {
        int r = m % n;
        if (r == 0) {
            return n;
        }
        m = n;
        n = r;
    }
}

#define min(a,b) (a<b?a:b)
#define abs(x,y) (x > y? x - y : y - x) 

uint32_t f(uint32_t y, uint32_t c, uint32_t n) {
    y = (y * y) % n;
    y += c;
    if (y < c)
        y += (std::numeric_limits<uint32_t>::max() - n) + 1;
    y %= n;
    return y;
}

//http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf
//pp:182-183
uint32_t rho(uint32_t n) {
    if (n == 1) return n;
    if (n % 2 == 0) return 2;

    uint32_t y = rand() % n; // y = x0
    uint32_t x;
    uint64_t ys = y;
    uint32_t c;
    do { c = rand() % n; } while (c == 0 || c == n - 2);
    uint32_t m = 1000;
    uint32_t d = 1;
    uint32_t q = 1;
    uint32_t r = 1;
    do
    {
        x = y;
        for (int i = 1; i <= r; i++)  y = f(y, c, n);
        int j = 0;
        do
        {
            ys = y;
            for (int i = 1; i <= min(m, r - j); i++)
            {
                y = f(y, c, n);
                q = q * abs(x, y) % n;
            }
            d = gcd(q, n);
            j += m;
        } while (j < r && d == 1);
        r *= 2;
    } while (d == 1);
    if (d == n)
    {
        do
        {
            ys = f(ys, c, n);//not:     ys = f(x,c,n);      
            d = gcd(abs(x, ys), n);
        } while (d == 1);
    }
    return d;
}