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 或 Tortoise and Hare 循环检测）。此时您将拥有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

Answer 2

编辑：
这对我有用，而不是被困在无限循环中：

#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
页：182-183
错误：ys = x * x % n;
正确：ys = ys * ys % n; // ys=f(ys)

请使用 stdint.h 中的 uint32_t 或 uint64_t ... 以获得良好而干净的编码风格：

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