浅谈 Miller-Robbin 与 Pollard Rho

前言

$Miller-Robbin$ 与 $Pollard Rho$ 虽然都是随机算法，不过用起来是真的爽。

$Miller Rabin$ 算法是一种高效的质数判断方法。虽然是一种不确定的质数判断法，但是在选择多种底数的情况下，正确率是可以接受的。

$Pollard Rho$ 是一个非常玄学的方式，用于在 $O(n^{1/4})$ 的期望时间复杂度内计算合数$n$的某个非平凡因子。

事实上算法导论给出的是 $O(sqrt p)$ ， $p$ 是 $n$ 的某个最小因子，满足 $p$ 与 $frac{n}{p}$ 互质。

但是这些都是期望，未必符合实际。但事实上 $Pollard Rho$ 算法在实际环境中运行的相当不错。

注：以上摘自洛谷。

Miller-Robbin

前置芝士

1. 费马小定理

• 内容：若 (varphi(p)=p-1,\,p>1)，则(a^{p}equiv apmod{p}) 或 (a^{p-1}equiv 1pmod{p},\,(a<p))

• 证明：戳这里

2. 二次探测定理

• 内容：如果 (varphi(p)=p-1,\,p>1,\,p>X) ，且(X^2equiv 1pmod{p})，那么(X=1 or p-1)

• 证明：

$ecause$ $X^2equiv 1pmod{p}$

$ herefore$ $p|X^2-1$

$ herefore$ $p|(X+1)(X-1)$

$ecause$ $p$是大于$X$的质数

$ herefore$ $p=X+1 or pequiv X-1pmod{p}$，即$X=1 or p-1$。

算法原理

由费马小定理，我们可以有一个大胆的想法：满足 $a^{p-1}equiv 1pmod{p}$ 的数字 $p$ 是一个质数。

可惜，这样的猜想是错误的，可以举出大量反例，如：$2^{340}equiv 1pmod{341}$，然鹅 $341=31*11$ 。

所以，我们可以取不同的 $a$ 多验证几次，不过，$forall a<561,\,a^{560}equiv 1pmod{561}$，然鹅 $561=51*11$ 。

这时，二次探测就有很大的用途了。结合费马小定理，正确率就相当高了。

这里推荐几个 $a_i$ 的值： $2，3，5，7，11，61，13，17$。用了这几个 $a_i$，就连那个被称为强伪素数的 $46856248255981$ 都能被除去。

主要步骤

• .将 (p-1) 提取出所有 (2) 的因数，假设有(s) 个。设剩下的部分为 (d)（这里提取所有(2)的因数，是为了下面应用二次探测定理） 。

• 枚举一个底数 (a_i) 并计算 (x=a_i^{d}pmod p)。

• 令 (y=x^{2}pmod p)，如果没有出现过 (p-1)，那么 (p) 未通过二次探测，不是质数。

• 否则，若底数已经足够，则跳出；否则回到第二步。

简易代码

```cpp #define ll long long ll p,a[]={2,3,5,7,11,61,13,17}; inline ll mul(ll a,ll b,ll mod) { ll ans=0; for(ll i=b;i;i>>=1) { if(i&1) ans=(ans+a)%p; a=(a<<1)%p; } return ans%p; } inline ll Pow(ll a,ll b,ll p) { ll ans=1; for(ll i=b;i;i>>=1) { if(i&1) ans=mul(ans,a,p); a=mul(a,a,p); } return ans%p; } bool Miller_Robbin(ll p) { if(p==2) return true; if(p==1 || !(p%2)) return false; ll d=p-1;int s=0; while(!(d&1)) d>>=1,++s; for(int i=0;i<=8 && a[i]Pollard Rho

大致流程

• 先判断当前数是否是素数（这里就可应用 (Miller-Robbin) ），如果是则直接返回

• 如果不是素数的话，试图找到当前数的一个因子（可以不是质因子）

• 递归对该因子和约去这个因子的另一个因子进行分解

如何找因子

一个一个试肯定是不行的。而这个算法的发明者采取了一种清奇的思路。（即采取随机化算法）

• 我们假设要找的因子为p

• 随机取一个 (x、y)，不断调整 (x) ，具体的办法通常是 (x=x*x+c)（c是随机的，也可以自己定）。

• 取 (p=gcd(y-x,n)) ，若(p in left(1,n ight)) ，则找到了一个因子，就返回。

• 如果出现 (x=y) 的循环，就说明出现了循环，并不断在这个环上生成以前生成过一次的数，所以我们必须写点东西来判环：我们可以用倍增的思想，让(y)记住(x)的位置，然后(x)再跑当前跑过次数的一倍的次数。这样不断让(y)记住(x)的位置，x再往下跑，因为倍增所以当(x)跑到(y)时，已经跑完一个圈。

• 同时最开始设定两个执行次数(i=1、k=2)（即倍增的时候用） ，每次取 (gcd) 时 (++i) ；如果 (i==k) ，则令 (y=x) ，并将 (k) 翻倍。

完整代码

• 题目：(Luogu4718 Pollard-Rho)算法 。

#include<cstdio>
#include<ctime>
#include<map>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('
')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
	TP f=1;x=0;register char c=getchar();
	for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
	for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
	x*=f;
}
template<typename TP>V print(TP x)
{
	if(x<0) x=-x,putchar('-');
	if(x>9) print(x/10);
	putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17,24251};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
template<typename TP>inline ll mul(TP a,TP b,TP p)
{
    ll ans=0;
    for(TP i=b;i;i>>=1)
    {
        if(i&1) ans=(ans+a)%p;
        a=(a<<1)%p;
    }
    return ans%p;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
	ll ans=1;
	for(TP i=b;i;i>>=1)
	{
		if(i&1) ans=mul(ans,a,p);
		a=mul(a,a,p);
	}
	return ans%p;
}
B Miller_Robbin(ll n)
{
	if(n<2) return false;
	if(n==2) return true;
	if(n%2==0) return false;
	ll d=n-1;int s=0;
	while(!(d&1)) d>>=1,++s;
	for(rg i=0;i<=8 && a[i]<n;++i)
	{
		ll x=Pow(a[i],d,n),y=0;
		F1(j,0,s-1)
		{
			y=mul(x,x,n);
			if(y==1 && x!=1 && x!=(n-1)) return false;
			x=y;
		}
		if(y!=1) return false;
	}
	return true;
}
inline ll Pollard_Rho(ll n)
{
	ll x,y,c,i,k;
	while(true)
	{
		ll x=rand()%(n-2)+1;
		ll y=rand()%(n-2)+1;
		ll c=rand()%(n-2)+1;
		i=0,k=1;
		while(++i)
		{
			x=(mul(x,x,n)+c)%n;
			if(x==y) break;
			ll d=Gcd(abs(y-x),n);
			if(d>1) return d;
			if(i==k) {y=x;k<<=1;}
		}
	}
}
V Find(ll n)
{
	if(n==1) return;
	if(Miller_Robbin(n)) {ans=max(ans,n);return;}
	ll p=n;
	while(n<=p) p=Pollard_Rho(p);
	Find(p);Find(n/p);
}
int main()
{
	read(T);srand(time(0));
	while(T--)
	{
		ans=0;
		read(n);Find(n);
		if(ans==n) puts("Prime");
		else print(ans),ed;
	}
	return 0;
}

**emmmm (cdots) **

这数据也太毒瘤了吧！！

看来要疯狂卡常了

优化1（不如叫做卡常？）

蛋定的分析一波，我们发现除了 $Pollard-Rho$ 是 $O(n^{1/4})$ 的期望时间复杂度外， $gcd$ 和龟速乘都是 $O(log N)$ 的。

虽然这种复杂度已经很优秀了，可对于本题的数据（(T≤350) 、 (1≤n≤10^{18})），还是太 (cdots)

所以我们要果断摒弃这种很 $low$ 的龟速乘，改用一种暴力溢出的快速乘：

• 简单原理： (a imes b  mod  p=a imes b−left lfloor frac{a imes b}{p} ight floor)

• 用 long double 来处理这个 (left lfloor frac{a imes b}{p} ight floor)

• 然后处理一下浮点误差就可以了。

• 模数较大时可能会出锅。

• 不过出锅概率很小 (cdots)

如下

```cpp inline ll mul(ll a,ll b,ll mod) { a%=mod,b%=mod; ll c=(long double)a*b/mod; ll ret=a*b-c*mod; if(ret<0) ret+=mod; else if(ret>=mod) ret-=mod; return ret; } ```

实践证明，战果辉煌：$6pts -> 94pts$ ！！！

优化2（正解）

虽然关于龟速乘的 $O(log n)$ 的恶劣影响得到了一定遏制，不过，我还是好想 $AC$ 啊！

通过办法1 ：打表 (cdots)

正确 $AC$ 姿势如下：

- 我们发现在 $Pollard-Rho$ 中如果长时间随机化而得不到结果，$gcd$带来的 $O(log n)$ 还是很伤肾的！！那有没有办法优化呢？答案是肯定的。

• 在生成(x)的操作中，龟速乘所模的数就是(n)，而要求的就是(n)的某一个约数，即现在的模数并不是一个质数。

• 根据取模的性质：如果模数和被模的数都含有一个公约数，那么这次模运算的结果必然也会是这个公约数的倍数。所以如果我们将若干个((y−x)) 相乘，因为模数是 (n) ，所以如果若干个((y−x))中有一个与(n)有公约数，最后的结果定然也会含有这个公约数。

• 所以可以多算几次((y−x))的乘积再来求(gcd) （一般连续算(127)次再求一次(gcd)）。

• 不过(cdots)，如果在不断尝试(x)的值时碰上一个环，就可能会还没算到(127)次就跳出这个环了，就无法得出答案；同时，可能(127)次计算之后，所有((y−x))的乘积都变成了(n)的倍数（即(prod_{i=1}^{127} {(y-x)} equiv 0 pmod{n}) ）

• 所以我们可以不仅在每计算(127)次之后求(gcd)、还要在倍增时（即判环时）求(gcd)，这样既维护了其正确性，又判了环！！

• 完整(AC)代码：

#include<cstdio>
#include<ctime>
#include<algorithm>
using namespace std;
#define rg register int
#define V inline void
#define I inline int
#define db double
#define B inline bool
#define ll long long
#define F1(i,a,b) for(rg i=a;i<=b;++i)
#define F3(i,a,b,c) for(rg i=a;i<=b;i+=c)
#define ed putchar('
')
#define bl putchar(' ')
template<typename TP>V read(TP &x)
{
	TP f=1;x=0;register char c=getchar();
	for(;c>'9'||c<'0';c=getchar()) if(c=='-') f=-1;
	for(;c>='0'&&c<='9';c=getchar()) x=(x<<3)+(x<<1)+(c^48);
	x*=f;
}
template<typename TP>V print(TP x)
{
	if(x<0) x=-x,putchar('-');
	if(x>9) print(x/10);
	putchar(x%10+'0');
}
int T;
ll n,ans;
ll a[]={2,3,5,7,11,61,13,17};
template<typename TP>inline ll Gcd(TP a,TP b) {return !b?a:Gcd(b,a%b);}
inline ll mul(ll a,ll b,ll mod)
{
	a%=mod,b%=mod;
	ll c=(long double)a*b/mod;
	ll ret=a*b-c*mod;
	if(ret<0) ret+=mod;
	else if(ret>=mod) ret-=mod;
	return ret;
}
template<typename TP>inline ll Pow(TP a,TP b,TP p)
{
	ll ans=1;
	for(TP i=b;i;i>>=1)
	{
		if(i&1) ans=mul(ans,a,p);
		a=mul(a,a,p);
	}
	return ans%p;
}
B Miller_Robbin(ll n)
{
	if(n<2) return false;
	if(n==2) return true;
	if(n%2==0) return false;
	ll d=n-1;int s=0;
	while(!(d&1)) d>>=1,++s;
	for(rg i=0;i<=8 && a[i]<n;++i)
	{
		ll x=Pow(a[i],d,n),y=0;
		F1(j,0,s-1)
		{
			y=mul(x,x,n);
			if(y==1 && x!=1 && x!=(n-1)) return false;
			x=y;
		}
		if(y!=1) return false;
	}
	return true;
}
inline ll Pollard_Rho(ll n)
{
	while(true)
	{
		ll x=rand()%(n-2)+1;
		ll y=rand()%(n-2)+1;
		ll c=rand()%(n-2)+1;
		ll i=0,k=1,b=1;
		while(++i)
		{
			x=(mul(x,x,n)+c)%n;
			b=mul(b,abs(y-x),n);
			if(x==y || !b) break;
			if(!(i%127) || i==k)
			{
				ll d=Gcd(b,n);
				if(d>1) return d;
				if(i==k) y=x,k<<=1;
			}
		}
	}
}
V Find(ll n)
{
	if(n<=ans) return;
	if(Miller_Robbin(n)) {ans=max(ans,n);return;}
	ll p=Pollard_Rho(n);
	while(n%p==0) n/=p;
	Find(p),Find(n);
}
int main()
{
	read(T);srand(time(0));
	while(T--)
	{
		ans=0;
		read(n);Find(n);
		if(ans==n) puts("Prime");
		else print(ans),ed;
	}
	return 0;
}