分解因子算法——Pollard $\rho$ 算法

攻击RSA密码体制最明显的方法就是分解大整数。关于大整数分解的算法有许多，常见的有试除法、$p\pm 1$算法、Pollard $\rho$算法、数域筛法、二次筛法等等。这篇分享的就是其中之一的Pollard $\rho$算法。

我们说因子分解，很多时候并不一定要把大整数$n$彻底分解成素数乘积，而是求出$n$的某个非平凡因子即可。求大整数的一个因子固然困难，但是Euclidean算法告诉我们求两个数的最大公因子是可以在$O(\log n)$的时间得到的。Pollard $\rho$算法利用的就是这个思想。

问题描述

设大整数为$n$，我们要找到$n$的一个因子$p$。

Pollard $\rho$ 算法

为了利用Euclidean算法找到$p$，我们设想存在两个整数$x，x^{\prime }\in {\mathbb{Z}}_n$，满足$x\equiv x^{\prime }(\mathrm{m}\mathrm{o}\mathrm{d}p)$，这样的话就有$p|(x-x^{\prime })$，又因为$p|n$，所以有$p\le gcd(x-x^{\prime },n)<n$——这句话的意思是：$gcd(x-x^{\prime },n)$一定是$n$的非平凡因子。这就完成了因子分解。问题归结为：如何找到满足条件的$x$和$x^{\prime }$。

一个朴素的想法就是：随机选择一个子集$X\subset {\mathbb{Z}}_n$，然后对所有不同的$x,x^{\prime }\in X$，判断$gcd(x-x^{\prime },n)$是否为1。这个方法能够成功当且仅当映射$x\to x(\mathrm{m}\mathrm{o}\mathrm{d}p)$在$X$中至少存在一个碰撞。根据生日攻击可以分析：如果$|X|\approx 1.71\sqrt{p}$时，至少存在一个碰撞的概率是50%，为了找到$X$中存在的这个碰撞，我们需要求$\left(\begin{array}{c}|X|\\ 2\end{array}\right)$次gcd。

Pollard $\rho$算法减少了求解gcd的次数。它并不是随机产生一个子集$X$，而是事先确定一个整系数多项式$f$，并取$x_1\in {\mathbb{Z}}_n$，令
$$x_{i+1}=f(x_i)\mathrm{m}\mathrm{o}\mathrm{d}n，i=1,2,3\cdots$$
这就产生了集合$X=\{x_1,x_2,\cdots ,x_m\}$。下面描述这个集合是如何减少gcd计算次数的。

核心思想：如果$x_i\equiv x_j(\mathrm{m}\mathrm{o}\mathrm{d}p)$，由于$f$是整系数多项式，所以$f(x_i)\equiv f(x_j)(\mathrm{m}\mathrm{o}\mathrm{d}p)$ 。又因为$p|n$，所以
$$x_{i+1}\mathrm{m}\mathrm{o}\mathrm{d}p=(f(x_i)\mathrm{m}\mathrm{o}\mathrm{d}n)\mathrm{m}\mathrm{o}\mathrm{d}p=f(x_i)\mathrm{m}\mathrm{o}\mathrm{d}p$$
同理有$x_{j+1}\mathrm{m}\mathrm{o}\mathrm{d}p=f(x_j)\mathrm{m}\mathrm{o}\mathrm{d}p$。因此$x_{i+1}\equiv x_{j+1}(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

进一步有：如果$x_i\equiv x_j(\mathrm{m}\mathrm{o}\mathrm{d}p)$，则对$\forall \delta \ge 0$ ,有$x_{i+\delta }\equiv x_{j+\delta }(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

记$l=j-i$，可知$x_{i^{\prime }}\equiv x_{j^{\prime }}(\mathrm{m}\mathrm{o}\mathrm{d}p)$如果$j^{\prime }\ge i^{\prime }\ge i$且$j^{\prime }-i^{\prime }\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}l)$。

这样序列X就可以看成带有一个“尾巴”
$$x_1\to x_2\to \cdots \to x_i(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
和一个“环”
$$x_i\to x_{i+1}\to \cdots \to x_j=x_i(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
看成图的话就像希腊字母$\rho$，这正是算法名字的由来。

对这个算法，我们需要计算$\left(\begin{array}{c}j-1\\ 2\end{array}\right)+i$次gcd才能得到第一个碰撞$x_i,x_j$。

进一步改进：我们只要求得一个碰撞即可，不一定是$x_i$和$x_j$，也可以是$x_{i+\delta }$和$x_{j+\delta }$。所以在找碰撞的时候，不必逐个逐个求gcd，可以跳跃着求。我们断言：在第一个“环”$x_i\to x_{i+1}\to \cdots \to x_j=x_i(\mathrm{m}\mathrm{o}\mathrm{d}p)$中，一定存在$i^{\prime }$满足$x_{2i^{\prime }}\equiv x_{i^{\prime }}(\mathrm{m}\mathrm{o}\mathrm{d}p)$。事实上，只要$2i^{\prime }-i^{\prime }=i^{\prime }=kl$，就有$x_{2i^{\prime }}\equiv x_{i^{\prime }}(\mathrm{m}\mathrm{o}\mathrm{d}p)$成立；而$i\le i^{\prime }=kl\le j-1$，所以最多找$j$次就可以得到$i^{\prime }$。

算法流程

Pollard $\rho$ (n,$x_1$)

$x=x_1，x^{\prime }=f(x)\mathrm{m}\mathrm{o}\mathrm{d}n$

$p=gcd(x-x^{\prime },n)$

while $p==1$

      //第i次迭代

​      $x=x_i,x^{\prime }=x_{2i}$

​       $x=f(x)\mathrm{m}\mathrm{o}\mathrm{d}p$

​       $x^{\prime }=f(x^{\prime })\mathrm{m}\mathrm{o}\mathrm{d}p$

​       $x^{\prime }=f(x^{\prime })\mathrm{m}\mathrm{o}\mathrm{d}p$

​       $p=gcd(x-x^{\prime },n)$

end

if p==n

​       return “failure”

else

​      return p

算法性能分析

假设$f=x^2+1$，则$j$的最大值为$\sqrt{p}$（因为$x_j$是第一个环中的数），所以最多需要$\sqrt{p}$次gcd计算。而$p<\sqrt{n}$，所以算法的期望复杂度是$O(n^{\frac{1}{4}})$。需要强调的是，这是一个启发式算法，不是严格的数学证明。算法可能会失败，此时是因为子集$X$中没有碰撞，需要重新选择初始值$x_1$或者选择不同的函数$f$产生新的集合$X$。

例子

设$n=7171$，$f(x)=x^2+1，x_1=1$。

计算集合$X$部分元素为：
$$1,2,5,26,677,6557,4105,6347,4903,2218,219,4936,4210,4560,4872,375,4377,4389,2016,5471,88\cdots$$
计算$gcd(x_1-x_2，n),gcd(x_2-x_4,n),gcd(x_3-x_6,n)\cdots$，发现$gcd(x_{11}-x_{22},n)=71$，这就找到了$n$的一个因子71。

参考书籍：Stinson D , 斯廷森, 冯登国. 密码学原理与实践[M]. 电子工业出版社, 2009.