C++

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

我改编我的代码的示例：

编辑

我已经按照 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

编辑：
这对我有用，而不是陷入无限循环：

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

旧：
根据改进的 monte carlo 分解算法：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;
}

Answer 2

它应该总是终止。当你到达算法中的那个点时，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] 范围内的数字。

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

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

C++ - Brent-Pollard rho algorithm infinite loop

algorithm

loops

factorization

编辑