Pollard's Rho Algorithm

参考资料

关于质因数分解与随机算法

我们假设要分解的数$n=pq$，其中$p\le q$，且$p$，$q$都为质数。这是最坏的情况。

肯定不能够$O(\sqrt n)$地去枚举$p$啦，这样的复杂度是无法承受的。我们考虑一下随机算法。我们随机一个小于$n$的数，它等于$p$（为$n$的某个约数）的概率为$1\over n$。那么，我们随机$n$次，期望会有$1$次，随机出来数是$p$。欸怎么好像还不如$O(\sqrt n)$……

那我们换个随机方式，随机$k$个数出来，把这$k$个数两两做差，看它们的差值是否为$p$（即是否为$n$的约数）。根据生日悖论，当$k$是$O(\sqrt n)$级别的时候，这$k$个数两两差值中，$p$的期望出现次数为$1$。但是把$O(\sqrt n)$个数拿来两两做差，复杂度是$O(n)$的，还是不如$O(\sqrt n)$的朴素算法。

“等于$p$”这个条件太严苛了，我们考虑放宽一点。假如随机出一个数$x$，我们只要$gcd(x,n)\not = 1$，就可以找到$n$的一个非平凡约数了（也就是这个$gcd$）。换言之，我们只需要得到一个$p$的倍数就可以了。那么根据生日悖论，当$k$达到$O(\sqrt p)$的时候，我们得到的这$k$个数中，期望有一对数，它们的差为$p$的倍数。但是这样检查的复杂度是$O(p)$的，对于$p=\sqrt n$（如果$n=10^{18}$，那么$p$将达到$10^9$）的情况仍然不够快，而且空间开销很大。

Pollard’s Rho algorithm

先看代码：

ll get(ll n,ll a)
{
	ll i=1,k=2,x=Rand()%n,y=x;// 这里不要取 Rand()%(n-1)+1 啥的，在n=4的时候会挂
	while(1)
	{
		++i,x=(Mul(x,x,n)+a)%n;
		ll d=gcd(Abs(y-x),n);
		if(d>1&&d<n) return d;
		if(x==y) return n;
		if(i==k) y=x,k<<=1;
	}
}

我们选择了一个伪随机数生成函数来生成一个随机数列，即$x_i = x_{i-1}^2 +a\mod n$。这里的$a$可以随意选取（不过算导说，$a=0,-2$应该被避免）。这个数列的取值，一定会陷入一个$\rho $形的环，即从某个位置开始，后面的元素的值是循环出现的。

我们发现，这个数列模$p$过后得到的数列，具有相同的递推式，即$x_i = x_{i-1}^2 + a\mod p$。那么相似地，这个模$p$意义下的序列也会陷入循环。现在我们的任务，是找到两个数$x_i,x_j$，使得它们在模$p$的意义下相等（也就是说，它们的差为$p$的倍数）并且它们在模$n$的意义下不相等。我们考虑在模$n$的这个数列的循环节上，固定一个位置，然后检查这个位置后面的$2^k$个位置，如果没有找到与它模$p$意义下相等的数——这个条件等价于没有找到一个数，它们的差与$n$的最大公约数不为$1$，也等价于没有找到一个数，它们在模$p$的那个数列的循环节上处于同一个位置——那么就让$k++$。如果在模$n$的这个数列的循环节上已经遍历完了一圈，这个时候仍然没有找到模$p$意义下与$y$相等的数，就直接返回，然后换一个种子或者换一个$a$重新计算。

这个算法的期望时间复杂度是$O(\sqrt {\sqrt n})$的，空间则是$O(1)$的。当然，前提是我们必须先判定，$n$不是素数，是可以分解的。这就需要Miller-Rabin素数测试啦。

上面只是找出一个约数的算法，要把一个大合数分解质因数，在外面套一个dfs就可以了。

// Pollard's rho algorithm
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll Abs(ll x){return x>0?x:-x;}
ll get(ll n,ll a)
{
	ll i=1,k=2,x=Rand()%n,y=x;
	while(1)
	{
		++i,x=(Mul(x,x,n)+a)%n;
		ll d=gcd(Abs(y-x),n);
		if(d>1&&d<n) return d;
		if(x==y) return n;
		if(i==k) y=x,k<<=1;
	}
}
vector<ll> fac;
void dfs(ll n)
{
	if(Isprime(n)){fac.PB(n); return;}
	ll p=n;
	while(p==n) p=get(n,Rand()%n);
	dfs(p),dfs(n/p);
}

这里有生成64位随机数的代码：

ull Rand()
{
	seed^=seed<<13,seed^=seed>>7,seed^=seed<<17;
	return seed;
}

模板：poj1811