3

下面是我对质数分解的 Pollard rho 算法的实现：

#include <vector>
#include <queue>

#include <gmpxx.h>

// Interface to the GMP random number functions.
gmp_randclass rng(gmp_randinit_default);

// Returns a divisor of N using Pollard's rho method.
mpz_class getDivisor(const mpz_class &N)
{
    mpz_class c = rng.get_z_range(N);
    mpz_class x = rng.get_z_range(N);
    mpz_class y = x;
    mpz_class d = 1;
    mpz_class z;

    while (d == 1) {
        x = (x*x + c) % N;
        y = (y*y + c) % N;
        y = (y*y + c) % N;
        z = x - y;
        mpz_gcd(d.get_mpz_t(), z.get_mpz_t(), N.get_mpz_t());
    }

    return d;
}

// Adds the prime factors of N to the given vector.
void factor(const mpz_class &N, std::vector<mpz_class> &factors)
{
    std::queue<mpz_class> to_factor;
    to_factor.push(N);

    while (!to_factor.empty()) {
        mpz_class n = to_factor.front();
        to_factor.pop();
        if (n == 1) {
            continue; // Trivial factor.
        } else if (mpz_probab_prime_p(n.get_mpz_t(), 5)) {
            // n is a prime.
            factors.push_back(n);
        } else {
            // n is a composite, so push its factors on the queue.
            mpz_class d = getDivisor(n);
            to_factor.push(d);
            to_factor.push(n/d);
        }
    }
}

它本质上是 Wikipedia 上伪代码的直接翻译，并且依赖 GMP 进行大数和素数测试。该实现运行良好，并且可以考虑素数，例如

1000036317378699858851366323 = 1000014599 * 1000003357 * 1000018361

但会窒息，例如

1000000000002322140000000048599822299 = 1000000000002301019 * 1000000000000021121

我的问题是：除了切换到更复杂的分解算法（例如Quadratic sieve ）之外，我能做些什么来改进这一点？

我知道一个改进可能是首先通过预先计算的素数进行一些试验除法，但这对于上述几个大素数的产品没有帮助。

我对有关改进基本 Pollard rho 方法的任何提示感兴趣，以使其处理仅包含几个主要因素的较大复合材料。当然，如果您在上面的代码中发现任何愚蠢之处，我也对它们感兴趣。

完整披露：这是一项家庭作业，因此一般提示和指针比完全编码的解决方案更好。通过这种非常简单的方法，我已经获得了作业的及格分数，但我当然想改进。

提前致谢。

4

1 回答 1

3

由于 Pollard，您正在使用 rho 算法的原始版本。Brent 的变体做了两处改进： Floyd 的龟兔赛跑算法被 Brent 开发的寻环算法所取代，并且 gcd 计算被延迟，因此它只在循环中每百次左右执行一次，而不是每次。但是这些变化只得到了很小的改进，可能是 25% 左右，并且不允许你考虑你正在谈论的大量数字。因此，您将需要一个更好的算法：SQUFOF 可能适用于您提到的大小的半素数，或者您可以实现二次筛或椭圆曲线方法。我在我的博客上讨论和实现了所有这些算法。

于 2013-11-02T15:49:10.987 回答