pollard’s rho算法(使用GMP库)

前言

整数分解是一个迷人的话题，至今为止也有不少的整数分解算法。有些算法是通用的整数分解算法，有些算法则是针对特殊的整数(比如之前介绍的pollard’s p-1算法)
下面我介绍一个通用的整数分解算法——pollard’s rho方法

注：这里的实现版本是我在密码学计算方法课程中学到的版本，和《算法导论》中的描述有一些区别，但算法的核心思想都是生日悖论和Floyd循环检测

算法描述

POLLARD_RHO
input: 合数n(n > 0)
output: n的一个因子

x = x₀(x₀为任意常数),
y = x² + 1
g = gcd(x – y, n)
While g == 1
x = x² + 1 mod n
y = (y²+1)²+1 mod n
End While

if 1 < g < n return g
else 失败(重新选择初始值x₀)

说明

pollard’s rho方法可以较快的找到合数的较小因子。这是一个比较简单的分解方法，它在分解素因子在10⁶附近的整数时，速度要比试除法快好几倍。但是它有一个不足，就是它是概率意义上的方法，必须对“坏情况”的出现有足够的警惕。除此之外，我们难以证明它的执行效果，从试验角度看，这可能认为是没有结果的。在比较重视正确性的应用中，这显然是一个致命的缺陷。

改进

在H.Cohen《A Course in Computational Algebraic Number Theory》一书中，作者认为应付“失效”状态的最好办法是使用一个不同于f(x) = x²+1的函数。f(y) = y²+c(c为常数)是一个合适的选择，尽管当c = 0, -2(mod n)时是“坏情况”。该书还断言，简单的改变x的初始值，也就是取不同x₀，不是聪明之举。

代码实现

//res为返回值，返回一个素因子
//n为待分解的合数
//p为x的初始值
void pollard_rho(mpz_t res, const mpz_t n, unsigned p)
{
    mpz_t x, y, g, t;
    mpz_init(x); mpz_init(y); mpz_init(g); mpz_init(t);
    mpz_set_ui(x, p);                                                   //x=p
    mpz_powm_ui(y, x, 2, n); mpz_add_ui(y, y, 1); mpz_mod(y, y, n);     //y=(x*x+1)%n
    mpz_sub(t, y, x);
    mpz_gcd(g, t, n);                                                   //g= gcd(y-x,n)

    while (mpz_cmp_ui(g, 1) == 0) {
        mpz_powm_ui(x, x, 2, n); mpz_add_ui(x, x, 1); mpz_mod(x, x, n); //x=(x*x+1)%n
        mpz_powm_ui(y, y, 2, n); mpz_add_ui(y, y, 1); mpz_mod(y, y, n); //y=(y*y+1)%n
        mpz_powm_ui(y, y, 2, n); mpz_add_ui(y, y, 1); mpz_mod(y, y, n); //y=(y*y+1)%n
        mpz_sub(t, y, x);
        mpz_gcd(g, t, n);                                               //g=gcd(y-x,n)

    }

    mpz_set(res, g);
    mpz_clear(x); mpz_clear(y); mpz_clear(g); mpz_clear(t);
}

运行结果

输入输出格式

输入
两行数据，每行一个整数
第一行为待分解的合数n
第二行为x的初始值p
输出
合数n的一个因子

样例输入
455459
5
样例输出
743

截图

附加代码

从该算法的适用范围来看，用gmp大数库显得有些多余，不过作为学习和熟悉gmp大数库的用法来说，这是值得的。下面的代码是我刚开始实现pollard’s rho 算法的原始版本(也是我提供给江智宇的版本)

int gcd(int a, int b)
{
    return (b == 0) ? a : gcd(b, a%b);
}

int pollard_rho(int n, int p)
{
    long long x = p;
    long long y = (x*x + 1) % n;
    int g = gcd(y-x+n, n);
    //printf("%15lld%15lld%15d\n", x, y, g);

    while (g == 1) {
        x = (x*x+1) % n;
        y = (y*y+1) % n;
        y = (y*y+1) % n;
        g = gcd(y-x+n, n);
        //printf("%15lld%15lld%15d\n", x, y, g);
    }

    return g;
}

总结

这个算法是我在密码学计算方法课程上学的其中一个整数分解方法，实现起来也非常简单，但是要弄清整个算法的来龙去脉，我估计我们学校没有学生可以做到。我对这个算法的了解程度也只是能够实现罢了。总而言之，算法分析是一个我可望而不可及的领域，过过瘾就可以了。
至于更加高级的因数分解算法——二次筛法(尽管核心思想很简单)，数域筛法，恐怕我这辈子都没能力写了。
围绕这个pollard’s rho算法还发生过有趣的事情。期末卷面考试有一题是用pollard’s rho算法分解455459，为了使得计算gcd的步骤最少，我们在考试前测试了许多初始值，最终发现当x=49或50时，只要计算两次gcd，一次循环(如下表)
n=455459 x=49

n=455459 x=50

最终我们选择了初始值为50作为试卷答案，因为124488这个数字比较顺口。
这是一个有趣的小伎俩。
不过话说回来，我对目前所学的数论算法都只是停留在了解层面，这些算法对于现在的我来说更像“魔法”，其背后的数学原理我并不特别清楚。