C++ - бесконечный цикл алгоритма Ро Брент-Полларда

Он всегда должен заканчиваться. К тому времени, когда вы дойдете до этой точки в алгоритме, x будет внутри цикла (поскольку y началось с x и прошло весь путь, чтобы обнаружить цикл, используя обнаружение цикла Брента). Значение ys началось с y, которое началось с x, поэтому оно также будет продолжаться по циклу и в конечном итоге снова встретится с x (см. Обнаружение цикла Флойда или Черепахи и Зайца). В этот момент у вас будет gcd(ys,x)==x, и этот последний цикл завершится.

В опубликованной реализации есть несколько ошибок, которые, как я полагаю, могут вызывать проблему:

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

Вы также можете заменить инициализацию c на

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

если вам нужно число в диапазоне [1,n-3].

Вот оригинал статьи Ричарда П. Брента: Улучшенный метод Монте-Карло Алгоритм факторизации

Редактировать:
это работает для меня и не попадает в бесконечный цикл:

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

Старый:
в соответствии с усовершенствованным алгоритмом факторизации Монте-Карло: http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf
pp:182-183
ошибка: ys = x * x % n;
правильно: ys = ys * ys % n; // ys=f(ys)

пожалуйста, используйте uint32_t или uint64_t... из stdint.h для хорошего и чистого стиля кодирования:

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