Miller-Rabin and Pollard-Rho

Miller-Rabin

她是一个素数判定的算法。

首先需要知道费马小定理

$$a^{p-1}\equiv1\pmod{p}\quad p\in prime$$

和二次探测定理

$$x=1 \ or \ p-1 \iff x^2\equiv1\pmod{p} \\ p\in prime,x\in[0,p)$$

证明

正确的运用平方差公式：

$$x^2\equiv1\implies (x+1)(x-1)\equiv 0\pmod{p}$$

然而对于素数 $p$，有（零因子）

$$ab\equiv0\implies a\equiv0\lor b\equiv0\pmod{p}$$

所以原命题得证。

Miller-Rabin（设要被检验的数为 $n$，$n\leqslant2 \ or \ 2\nmid n$ 的情况请特判）:

1. 将 $n-1$ 拆成 $2^{cnt}\times m\quad (2\nmid m)$ 的形式，并随机选出一个整数 $a$ ，构造此 Miller 序列：

$$\begin{aligned} M_i&=a^{2^i\times m}\pmod{n}\quad (0\leqslant i\leqslant cnt) \\ &=\begin{cases} a^m\pmod{n} & i=0 \\ M_{i-1}^2\pmod{n} & i>0 \end{cases} \end{aligned}$$

2. 由于费马小定理，若 $M_{cnt}\ne1$，则 $n$ 一定不是素数；若 $M_{cnt}=1$，则 $n$ 可能是素数。

3. 若 $n$ 没有被第二步筛掉，即还有机会是素数，则研究 $M_{cnt-1}$：若
   $M_{cnt-1}=1 \ or \ n-1$，则二次探测定理成立， $n$ 可能是素数；否则二次探测定理不成立， $n$ 一定不是素数。

4. 若 $M_{cnt-1}=n-1$ ，则接下来不能用二次探测定理判定了，直接返回 $n$ 为素数。

5. 若 $n$ 没有被第三步筛掉，即还有机会是素数，则研究 $M_{cnt-2}$：若
   $M_{cnt-2}=1 \ or \ n-1$，则二次探测定理成立， $n$ 可能是素数；否则二次探测定理不成立， $n$ 一定不是素数。

6. 递归地使用第四第五两步，若到了 $M_0$ 还没有出结果，即 $M_i=1\quad (0\leqslant i\leqslant cnt)$ ，则返回 $n$ 为素数。

7. 重复多次上述过程，即 $a$ 多次随机。只有当每次结果都返回是素数时，最终判定结果才是素数，否则不是素数。

举个例子,

若要判定 $561=3\times11\times17$ 是否为素数:

$561-1=560=2^4\times35,cnt=4,m=35$。$a$ 取 $8$，则构造的 Miller 序列分别为：

$$M_0=a^{2^0\times 35}\pmod{561}=461$$

$$M_1=a^{2^1\times 35}\pmod{561}=463$$

$$M_2=a^{2^2\times 35}\pmod{561}=67$$

$$M_3=a^{2^3\times 35}\pmod{561}=1$$

$$M_4=a^{2^4\times 35}\pmod{561}=1$$

发现 $M_2$ 和 $M_3$ 不满足二次探测定理，断定了 $561$ 不是素数，事实确实如此。

点击查看代码

inline bool MR(ll n){//Miller-Rabin
	#define Times 80//检验次数 
	if(n<2 || !(n&1))return (n==2);//特判
	ll a,b,m=n-1,cnt=0;//cnt:2因子的个数 
	while(!(m&1)){
		cnt++;
		m>>=1;
	}
	uid(R,1,n-1);//随机数init,这里用mt19937
	For(i,1,Times){
		a=Rand(R);//随机值
		a=pw(a,m,n);//快速幂
		For(j,1,cnt){
			b=a*a%n;//此时b为M_j,a为M_{j-1}
			if(b==1 && a!=1 && a!=n-1)//二次探测定理都不行qwq 
				return false;
			a=b;
		}
		if(b!=1)//费马小定理都不行qwq 
			return false;
	}
	return true; 
}

发现这样常数太大，若将 Times 调小，则有些素数会误判。

优化

选取几个特定的优秀的素数作为 $a$，此时还要特判这些素数

点击查看代码

inline bool mr(ll n,ll c){//用c来筛n 
    ll k=n-1,cnt=0;
    while(!(k&1)){
    	k>>=1;
    	cnt++;
	}
	ll x=pw(c,k,n),y=x;
	if(x==1) return true;
    while(cnt--){
    	y=pw(x,2,n);
    	if(y==1 && x!=n-1 && x!=1) return false;
		x=y;
    	if(x==1) return true;
    }
    return false;
}
vector<ll> pr={2,3,5,7,13331,998244353}; //素数选得好，AC不能少 
inline bool MR(ll x){
    if(x<2) return false;
    for(ll i:pr) if(x==i) return true;
    for(ll i:pr) if(!mr(x,i)) return false;
    return true;
}

Pollard-Rho

她是一个分解大合数的算法。

如果是你，你会如何找到某个大整数的一个非1、非自身的因子？

算法1：试除法

从 $2$ 到 $\sqrt{n}$ 枚举，判断是否为 $n$ 的因子。

时间复杂度 $O(\sqrt{n})$

算法2：随机数法

从 $2$ 到 $n-1$ 枚举，判断是否为 $n$ 的因子。

这个算法看上去像是来搞笑的，简直和猴子排序不谋而合，然而它正是Pollard-Rho算法的基础。

最坏情况下的期望时间复杂度 $O(n)$

算法3：2的优化

从判断是否为 $n$ 的因子改成判断是否与 $n$ 不互素，之后取 $\gcd$ 即可。

最坏情况为：

$$n=p^2\quad p\in prime$$

此时期望时间复杂度 $O(\sqrt{n}\log n)$ （$\log$ 是 $\gcd$ 带来的）

为了更好得介绍下面的内容，先引入一堆——

生日悖论（其实不是悖论）

这只是有悖直觉而已 OvO

一个房间里有23个人，则他们中有两人生日相同的概率超过一半（不考虑闰年）。

证明从略。

生日悖论启示我们，如果我们不断在某个范围内生成随机整数，很快便会生成到重复的数，期望大约在根号级别。

对于一个 $[1,N]$ 内整数的随机生成器，生成序列中第一个重复数字的位置的期望为

$$\sqrt{\dfrac{\pi N}{2}}$$

但这件事意义并没有那么大。正如虽然生日悖论是正确的，但你不一定能在班上遇到和自己生日相同的人，因为这个高概率是在两两比较下才成立的。对这些数两两进行验证，复杂度立刻退化，并没有什么进步。所以我们需要一些技巧。

伪随机数序列

设 $f(x)=x^2+c$（可能受曼德勃罗集的启发吧）

则构造序列

$$P(i)=\begin{cases} 0 & i=0 \\ f(P(i-1)) & i>0 \end{cases}$$

和

$$P(i)_n=P(i)\pmod{x}$$

由于 $P(i)_n$ 的取值只由上一个决定（马尔可夫链）（$n,c$ 看作常数）且值域有限，则此序列一定会形成环。

举个例子：

取 $n=1001=7\times11\times13,c=2$ 时，序列 $P(i)_n$ 的走向为：

生成的序列常常形成这样的 $\rho$ 形，这也是为什么 Pollard 把这个算法命名为 $\rho$ (rho) 算法。

算法4：Pollard-Rho

转换一下原问题至：不断求出 $n$ 的一个非 $1$ 因数，直到此因数不为 $n$ 即可。

接下来考虑如何求出 $n$ 的一个非 $1$ 因数。

其实就是找 $P$ 序列中两数 $x,y$（不一定相邻)，使得 $\gcd(|x-y|,n)>1$

我们称整个为 $\rho$，环（此图中的 $445\to830\to214\to753\to445$）为 $o$，$o$ 的环长为 $len(o)$

这伪随机数序列有什么好处呢？其实，这个伪随机数生成器生成的数具有一个性质：

若 $\rho$ 上有两个数 $i,j$ 位置相差 $d\quad(len(o) \nmid d)$ （$e.g.$ 图中 $38$ 与 $214$ 相差 $3$；$445$ 与 $753$ 相差 $1$ 或 $3$ 均可），且

$$\gcd(|i-j|,n)>1$$

则在 $o$ 上 $\forall x,y$ 位置相差 $d$ ,有

$$\gcd(|x-y|,n)>1$$

证明

$len(o) \nmid d \implies \rho$ 上位置相差 $d$ 的数不可能相同

$$\gcd(|i-j|,n)>1$$

$$\implies \gcd(|i-j|\times|i+j|,n)>1$$

$$\implies \gcd(|i-j|\times|i+j|,n)>1$$

$$\implies \gcd(|i^2-j^2|,n)>1$$

$$\implies \gcd(|f(i)-f(j)|,n)>1$$

$$\implies \dots$$

（递推）

$$\implies \gcd(|x-y|,n)>1$$

换言之，任意选取 $o$ 上距离为 $i\quad(0<i<len(o))$ 的两点 $x_i,y_i$，判断 $\gcd(|x_i-y_i|,n)$ 是否为 $n$ 的非 $1$ 因数。

这些二元组基本上（因为距离为 $len(o)$ d 倍数的没算，由于这是概率算法，这些可忽略）代表了 $\rho$ 上的所有二元组。

所以我们在一个 $o$ 上，对于每种距离 $i$，任意找一组 $x_i,y_i$ 判断是否有非 $1$ 因子即可。

若没找到，则 $f(x)$ 中 $c$ 换一个随机值构造 $\rho$ 重复上述步骤即可。

那怎么快速枚举环上每种距离的两点呢？——

Floyd判环算法（龟兔赛跑算法）

形象理解一下，就是把 $\rho$ 看作飞行棋棋盘，有两只小灰鸡从 $0$ 开始，慢的（称为 $turtle$）每次走一格，快的（称为 $rabbit$）每次走两格。

即初始

$$turtle=rabbit=0$$

每次

$$turtle:=f(turtle),rabbit:=f(f(rabbit))$$

且判断 $\gcd(|rabbit-turtle|,n)$

由于 $rabbit$ 与 $turtle$ 之间的距离先逐渐增大，之后在环上相遇，所以两者之间的距离几乎可以取到所有可能的值。

期望复杂度 $O(n^{\frac{1}{4}}\log n)$

点击查看代码

inline ll f(ll x,ll c,ll n){
    return (/*__int128*/(lll)x*x+c)%n;
}
inline ll PR(ll n,ll c){
    ll t=f(0,c,n);//turtle 
	ll r=f(f(0,c,n),c,n);//rabbit
    while(t!=r){
        ll d=gcd(abs(t-r),n);
        if(d>1)
            return d;
        t=f(t,c,n);
		r=f(f(r,c,n),c,n);
    }
    return n;//没有找到,重新调整参数c
}

据说这个方法分解不了 $4$，所以最好特判一下。

时间复杂度的证明

设 $d|n\quad(1<d\leqslant\sqrt{n})$

因为 $P(i)_n$ 和 $P(i)_d$ 可近似看为随机序列，根据生日悖论可以推出其出现循环的期望步数分别为

$$\sqrt{\dfrac{\pi n}{2}},\sqrt{\dfrac{\pi d}{2}}$$

由于 $n>d$ ，所以有极大的概率当第一次 $P(rabbit)_d=P(turtle)_d$ 时，$P(rabbit)_n\ne P(turtle)_n$（即 $P_d$ 比 $P_n$ 更早进入循环)

此时

$$\begin{aligned} |P(rabbit)_n-P(turtle)_n|&=|(k_{rabbit}d+P(rabbit)_d)-(k_{turtle}d+P(turtle)_d)| \\ &=|(k_{rabbit}-k_{turtle})d+(P(rabbit)_d-P(turtle)_d))| \\ &=|(k_{rabbit}-k_{turtle})d|\end{aligned}$$

为 $d$ 的倍数，于是求 $\gcd$ 即可求出 $d$ 这个因子。

因为 $P_d$ 的循环期望步数 $= \sqrt{\dfrac{\pi d}{2}} \approx \sqrt{d}$，而 $d\leqslant\sqrt{n}$，再乘上 $\gcd$ 的时间复杂度 $O(\log n)$，则最后是 $O(n^{\frac{1}{4}}\log n)$

优化：倍增距离

这个方法复杂度很低了，但是这个 $\log$ 还是看得人不爽，我们考虑如何去掉它。

我们想到

$$\gcd(x,n)>1 \implies \gcd(kx,n)>1 \quad (k\in \mathbb{N^+})$$

所以我们可以减少求公因数的次数，即先把一些待选数乘起来，再统一与 $n$ 求公因数。

我们可以每隔 $1,2,4,8,16,\dots$ 个数（依次增加）求一次公因数，这样只需要求期望 $\log(n^{\frac{1}{4}})$ 次公因数。

具体做法改为（伪代码）：

1. $turtle:=rabbit:=k:=0$

2. $mul:=1$

3. $\operatorname{For}(i:1\to 2^k) rabbit:=f(rabbit),mul:=mul \times |rabbit-turtle|$

4. $\operatorname{check\_gcd}(mul,n)$

5. $turtle:=rabbit$

6. $k:=k+1$

7. $\operatorname{goto}(step \ 2)$

即 $turtle$ 开启了瞬移到 $rabbit$ 身边的功能，冷却时间倍增。

总期望时间复杂度为 $O(n^{\frac{1}{4}}+\log(n^{\frac{1}{4}})\log(n))\approx O(n^{\frac{1}{4}})$

优化：固定步数

倍增距离在常数上有个缺点，就是到后面间隔很大，可能已达成目标却迟迟无法退出。

所以我们也可以龟兔赛跑算法（不带瞬移）每隔固定步数 $C\approx \log n$ （取 $128$ 较优）就求一次公因数。

注意：这里的固定步数不是指固定 $turtle,rabbit$ 之间的距离，而是与Floyd判环算法/龟兔赛跑算法相同，只不过将相邻的 $\gcd$ 询问堆在一起处理。

总期望时间复杂度为 $O(n^{\frac{1}{4}}+\dfrac{n^{\frac{1}{4}}\log n}{C})\approx O(n^{\frac{1}{4}})$

网上很多模板会把这两种方法结合，即倍增取距离，但又规定一个上限。本人本地测试了一下，似乎跟固定步数的时间差别不大，但 Luogu 上好像只有结合的方法能过，固定步数TLE调不出来 qwq。

点击查看代码

inline ll gcd(ll x,ll y){//定义 gcd(0,x)=x 
	if(x<y) swap(x,y);
	if(y==0) return x;
	return gcd(y,x%y);
}
inline ll f(ll x,ll c,ll n){
    return ((lll)x*x+c)%n;
}
inline ll PR(ll n,ll c){
    ll t;//turtle
	ll r=0;//rabbit 
    ll val;//累乘 
    ll d;//因数 
    for(int len=1;;len<<=1){//len:这一阶段的长度 （倍增） 
    	t=r;
    	val=1;
        For(cnt,1,len){//rabbit跑 
            r=f(r,c,n);
            val=(lll)val*abs(r-t)%n;
            if(cnt%127==0){
                d=gcd(val,n);
                if(d>1)
                    return d;
            }
        }
        d=gcd(val,n);
        if(d>1)
            return d;
    }
}

求最大因数（ Luogu P4718 ）

既然可以求一个因数，自然就可以求最大因数，这需要一个递归即可。

点击查看代码

ll ans;
void work(ll n){//求n的最大素因数->ans
	if(n<=ans || n<2) return ;
	if(MR(n)){//特判素数 
		ckmx(ans,n);
		return ;
	}
	ll d=n;
	uid(R,1,n-1);
	while(d==n) d=PR(n,Rand(R));//找到一个因子（不一定素） 
	while(n%d==0) n/=d;
	work(n);
	work(d);
}

完整呆码