C++ - Brent-Pollard rho 算法无限循环

Question

我有以下功能。我从两个站点获得它并尝试将其改编成我自己的，但效果不佳。

当我测试 unsigned long int max - 2 或 to 但它作为一个数字 4294967293 时，它会将以下代码放入无限循环中，其中 ys 将继续返回相同的值。

我对因式分解算法非常陌生，如果能逐步帮助我了解为什么会出现无限循环，那就太好了。

以下代码只是我的“rho”函数。我有另一个名为 gdc 的函数，它与所有其他 gdc 递归函数相同。

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

我改编我的代码的示例:

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

编辑

我已经按照 Amd 的建议完成并整理了我的代码，并将重复的行移到了辅助函数中。但是，对于底部附近的 d==n 部分，我仍然遇到无限循环。出于某种原因，f(ys) 最终基本上返回了它之前返回的相同内容，因此它不断循环通过一系列值。

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

它应该总是终止。当你到达算法中的那个点时，x 将在循环内(因为 y 从 x 开始并一直进行周围使用布伦特循环检测来检测循环)。值 ys 从 y 开始，从 x 开始，因此它也会继续循环并最终与 x< 相遇 再次(参见 Floyd 或龟兔赛跑检测)。此时您将得到 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] 范围内的数字。

这是 Richard P. Brent 的原始论文:An Improved Monte-Carlo Factorization Algorithm