Miller Rabin & Pollard Rho & Pollard P-1 & SQUFOF

Miller Rabin

费马素性检验

费马小定理，当 \(p\) 是质数时有 \(\forall a\in[1,p-1],a^{p-1}\equiv 1\pmod p\)。

因此可以在 \([0,p-1]\) 选几个 \(a\) 检验 \(a^{p-1}\equiv1\pmod p\)。

但是存在一些数（卡迈克尔数） 也满足 \(\forall a\in[1,p-1],a^{p-1}\equiv 1\pmod p\)。

---

二次探测

二次探测定理，当 \(p\) 是奇质数时有 \(x^2\equiv1\pmod p\) 的解只有 \(x\equiv\pm 1\pmod p\)。

---

Miller Rabin

设选取 \(a\) 作底数检验 \(p\)。

令指数 \(d=p-1\)。

• 若 \(a^{p-1}\not\equiv 1\pmod p\)，则 \(p\) 为合数。
• 有 \(a^d\equiv 1\pmod p\)。
• 令 \(d\gets\dfrac d 2\)，因为 \((a^{\frac d 2})^2\equiv 1\pmod p\)，所以应有 \(a^{\frac d 2}\equiv\pm 1\pmod p\)。
• 若 \(a^d\equiv -1\pmod p\)，退出（可能是质数）；否则重复前一步。

底数表：

\(2^{32}\) 内：\(\{2,7,61\}\)。

\(2^{64}\) 内：\(\{2,325,9375,28178,450775,9780504,179265022\}\)。

---

实现

\(O(\log^2 p)\)：

#define mul(x,y,m) ((__int128)(x)*(y)%m)
inline long long qpow(long long x,long long y,long long mod){
	long long res=1;
	while(y)
	{
		if(y&1){
			res=mul(res,x,mod);
		}
		y>>=1;
		x=mul(x,x,mod);
	}
	return res;
}
inline bool MR(int a,long long p){
	long long d=p-1;
	if(a%p==0){
		return true;
	}
	if(qpow(a,d,p)^1){
		return false;
	}
	while(~d&1)
	{
		d>>=1;
		long long cur=qpow(a,d,p);
		if(cur==p-1){
			return true;
		}
		if(cur^1){
			return false;
		}
	}
	return true;
}
inline bool Miller_Rabin(long long p){
	if(p<=1){
		return false;
	}
	static int tb[]={2,325,9375,28178,450775,9780504,179265022};
	for(long long a:tb)
	{
		if(!MR(a,p)){
			return false;
		}
	}
	return true;
}

\(O(\log p)\)：

inline bool MR(int a,long long p){
	long long d=p-1;
	if(a%p==0){
		return true;
	}
	int c=0;
	while(~d&1)
	{
		c++;
		d>>=1;
	}
	long long cur=qpow(a,d,p);
	if(cur==1||cur==p-1){
		return true;
	}
	for(int i=1;i<=c;i++)
	{
		cur=mul(cur,cur,p);
		if(cur==1){
			return false;
		}
		if(cur==p-1&&i^c){
			return true;
		}
	}
	return (cur==1);
}

---

Pollard Rho

生日悖论

假设有理想的 \([1,n]\) 的随机数生成器，期望生成 \(\sqrt\dfrac{\pi n}{2}\approx O(\sqrt n)\) 次生成出两个相同的数。

设待分解的合数为 \(n\)，其最小质因子为 \(p\le \sqrt n\)，那么期望 \(O(\sqrt p)=O(n^{\frac 1 4})\) 次生成出两个数 \(\bmod\ p\) 同余。

---

Pollard Rho

考虑用一种优美的伪随机方式，生成数列 \(\{a\}\)，\(a_i=({a_{i-1}}^2+c)\bmod n\)，即用 \(f(x)=x^2+c\) 生成。

期望 \(O(n^{\frac 1 4})\) 次出现 \(a_i\equiv a_j\pmod p\)，那么 \({a_{i-1}}^2\equiv{a_{j-1}}^2\pmod p\)，即 \((a_{i-1}-a_{j-1})(a_{i-1}+a_{j-1})\equiv 0\pmod p\)。

虽然这个性质似乎没什么用。

现在需要找到 \(f(x)\bmod p\) 的环，判定方式为 \(\gcd(|a_i-a_j|,n)>1\)。

注意有可能全程判不出来，这时候换一个 \(c\) 重新来即可。（有些 corner case 所有 \(c\) 都不行，提前除掉 \(2\) 即可）

---

Floyd 判环

设两个数 \(s=a_0,t=a_1\)，每次 \(s\gets f(s),t\gets f(f(t))\)，用 \(s,t\) 判环（\(\gcd(|s-t|,n)>1\)）。

只需要 \(O(len)\) 次判断，\(len\) 为环长。

单次有 \(\gcd\) 的 \(\log n\)，时间复杂度 \(O(n^{\frac 1 4}\log n)\)。

#include<random>
mt19937_64 rnd(372409);
inline long long gcd(long long a,long long b){
	long long az=__builtin_ctzll(a),bz=__builtin_ctzll(b),z=(az>bz)?bz:az,t;
	b>>=bz;
	while(a)
	{
		a>>=az;
		t=a-b;
		az=__builtin_ctzll(t);
		b=a<b?a:b;
		a=t<0?-t:t;
	}
	return b<<z;
}
inline long long Pollard_Rho(long long n){
	long long c=rnd()%n+1;
	#define f(x) (((__int128)(x)*(x)+c)%n)
	long long s=f(0),t=f(s);
	while(s^t)
	{
		long long g=gcd(abs(s-t),n);
		if(g>1){
			return g;
		}
		s=f(s);
		t=f(f(t));
	}
	return n;
}

---

Brent 判环

Floyd 判环的常数优化。

设 \(k=1\)，每轮 \(t\) 向前 \(2^k\) 次判环，判不出则 \(s\gets t\)，\(k\gets k+1\)，进入下一轮。

inline long long Pollard_Rho(long long n){
	long long c=rnd()%(n-1)+1;
	#define f(x) (((__int128)(x)*(x)+c)%n)
	long long s,t=f(0);
	for(int k=1;;k<<=1)
	{
		s=t;
		for(int i=0;i<k;i++)
		{
			t=f(t);
			long long g=gcd(abs(s-t),n);
			if(g>1){
				return g;
			}
		}
	}
	return n;
}

话说为啥我测出来比 Floyd 判环慢啊？

---

倍增优化

若 \(\gcd(x,n)>1\)，那么 \(\gcd(xy,n)>1\)。

记录 \(|s-t|\) 的乘积 \(v\)，每 \(O(\log n)\) 个判一下 \(\gcd(v,n)>1\)。

复杂度就能优化到近似 \(O(n^{\frac 1 4})\)。

inline long long Pollard_Rho(long long n){
	long long c=rnd()%(n-1)+1;
	#define f(x) (((__int128)(x)*(x)+c)%n)
	long long s,t=0;
	for(int k=1;;k<<=1)
	{
		s=t;
		long long v=1;
		for(int i=1;i<=k;i++)
		{
			t=f(t);
			v=mul(v,abs(s-t),n);
			if(!(i&0x7f)){
				long long g=gcd(v,n);
				if(g>1){
					return g;
				}
			}
		}
		long long g=gcd(v,n);
		if(g>1){
			return g;
		}
	}
	return n;
}

---

Pollard P-1

设 \(n-1\) 的最大质因子是 \(B\)，那么 \(\gcd(2^{B!}-1,n)\) 是 \(n\) 的一个因数。

费马小定理，\(a^{t(p-1)}-1=kp\)，也就是 \(p-1|t\rArr p|a^{t}-1\)。

设 \(n=pq\)，若存在底数 \(B\) 使得 \(p|a^B-1\) 且 \(q\nmid a^B-1\)，那么 \(\gcd(2^{B!}-1,n)\) 分解 \(n\)。

\(B\) 就是 \(n-1\) 的最大质因子，找到的就是 \(p,q\) 中 \(B-\mathrm{smooth}\) 的那个。

inline long long Pollard_P(long long n){
	long long sqr=(long long)(sqrtl(n)+1e-6);
	if(sqr*sqr==n){
		return sqr;
	}
	long long a=2;
	for(int j=2;;j++)
	{
		a=qpow(a,j,n);
		long long g=gcd(a+n-1,n);
		if(g>1){
			return g;
		}
	}
	return n;
}

现在还没找到复杂度分析，反正超级慢就是了。

---

SQUFOF

inline long long SQUFOF(long long n){
	if(~n&1){
		return 2;
	}
	long long sqr=(long long)(sqrtl(n)+1e-6);
	if(sqr*sqr==n){
		return sqr;
	}
	const long long tb[]={1,3,5};
	for(int k:tb)
	{
		if(k^1&&!(n%k)){
			return k;
		}
		long long K=k*n,R=(long long)(sqrtl(K)+1e-6);
		long long b=0,P=R,Q=1;
		#define trans(){\
			swap(b,Q);\
			long long tmp=P,q=(R+P)/b;\
			P=q*b-P;\
			Q+=q*(tmp-P);\
		}
		long long lim=(long long)(sqrtl(R<<1)+1e-6);
		lim<<=2;
		trans();
		Q=(K-P*P)/b;
		for(int i=2;i<lim;i+=2)
		{
			trans();
			sqr=(long long)(sqrtl(Q)+1e-6);
			if(sqr*sqr==Q){
				long long _b=b,_P=P,_Q=Q;
				P=-P;
				Q=sqr;
				trans();
				Q=(K-P*P)/b;
				do{
					sqr=P;
					trans();
				}
				while(P^sqr);
				sqr=gcd(n,sqr);
				if(sqr>1){
					return sqr;
				}
				b=_b;
				P=_P;
				Q=_Q;
			}
			trans();
		}
	}
	return n;
}