Pollard Rho 和 Miller Rabin

Miller Rabin

适用问题

判断一个 \(2^{100}\) 以内的数是否是质数。

算法思想

选取若干个数,依次以它们作为底数进行费马探测和二次探测,如果都说是质数,它就是质数。
底数的选择:前 12 个质数(适用于 \(2^{78}\) 以内的质数判断)。

费马探测

\({a^{p-1}\not\equiv 1\pmod p}\),则 \(p\) 不是质数。

二次探测

\(d\) 是偶数,若 \({a^d\equiv 1\pmod p}\not\Rightarrow {a^{\frac d2}\equiv ±1}\),则 \(p\) 不是质数。

可以先令 \(d=p-1\),迭代 \(d\gets d/2\) 直到 \(d\) 为奇数或 \(a^d\equiv -1\pmod p\),中途若遇到不是 \(±1\) 的就返回假。

代码

bool check(int p,int b){
	int d=p-1,tmp=qp(b,d,p);
	if(tmp!=1)return 0;
	while(~d&1){
		if(d>>=1,(tmp=qp(b,d,p))==p-1)return 1;
		else if(tmp!=1)return 0;
	}
	return 1;
}
bool miller(int p){
	if(p<40){//essential
		for(int i=0;i<12;i++)if(pr[i]==p)return 1;
		return 0;
	}
	for(int i=0;i<12;i++)if(!check(p,pr[i]))return 0;
	return 1;
}

Pollard Rho

适用问题

对一个 \(2^{100}\) 以内的数进行质因数分解。

算法思想

每次找到 \(n\) 的任意一个非平凡(\(\notin\{1,n\}\))因子 \(p\),并当 \(p|n\) 时迭代 \(n\gets n/p\),然后递归分解 \(p\) 和现在的 \(n\),以获得有关质因数分解的信息,直到用 Miller Rabin 判断 \(n\) 为质数时回溯。

如何求出任意一个非平凡因子:
生成伪随机数列 \(a=\{0,c,c^2+c,(c^2+c)^2+c,((c^2+c)^2+c)^2+c,...\pmod p\}\),由于是伪随机数列,不可避免地会在某一时刻进入环,即 \(a_i=a_j(i>j)\)。选择若干组 \((a_A,a_B)\),判断是否 \(\gcd(|a_A-a_B|,p)>1\),若是,则 \(\gcd(|a_A-a_B|,p)\) 符合条件。根据生日悖论,选取数对判断的成功率远远高于选取独立元素 \(\gcd(a_A,p)\) 判断的成功率。

写法

Floyd 判环

取两个数 \(x_0,y_0\),记伪随机函数为 \(f(t)=(t^2+c)\bmod p\),每次迭代令 \(x_0=f(x_0),y_0=f(f(y_0))\),判断 \((x_0,y_0)\) 数对。
可以证明此算法可在可观的时间内退出。

倍增法

假如走若干轮,第 \(k\) 轮走 \(2^k\) 步,这一轮开始时的数是 \(s\),令 \(t=s\),并开始迭代(走步)\(t=f(t)\),每走一步判断 \((s,t)\) 数对。
调用 __gcd 函数的时间复杂度为 \(\log n\),考虑减少调用。
考察 \(\gcd(|s-t_1|,p)>1\),则 \(\gcd(|s-t_1t_2|,p)>1\),进一步地,\(\gcd(|s-\prod_{j=1}^rt_j\bmod p|,p)>1\),故考虑过一段定长时间判一次。取阈值为 \(127\)
倍增法的效率高于 Floyd 判环。

代码(倍增法)

int pollard(int x){
	int c=rand()%(x-1)+1,s=0,t=0,d;
	for(int goal=1;;goal<<=1,s=t){
		int prod=1;
		for(int i=1;i<=goal;i++){
			t=((lll)t*t+c)%x;
			prod=(lll)prod*abs(t-s)%x;
			if(i%127==0){if((d=__gcd(prod,x))>1)return d;}
		}
		if((d=__gcd(prod,x))>1)return d;
	}
	return x;
}
void fac(int n){
	if(n<2)return;//还可加入其他剪枝
	if(miller(n)){Do something;return;}
	int p=0,k=0;
	while((p=pollard(n))==n);
	while(n%p==0)n/=p,k++;
	fac(p),fac(n);
}

习题:[POI2010]Divine Divisor

#include <bits/stdc++.h>
#define int long long
using namespace std;
typedef __int128 lll;
inline int read(){
	register char ch=getchar();register int x=0;
	while(ch<'0'||ch>'9')ch=getchar();
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
	return x;
}
inline int qp(int a,int b,int p){
	int c=1;for(;b;b>>=1,a=(lll)a*a%p)if(b&1)c=(lll)c*a%p;
	return c;
}
int n,max_factor,pr[13]={2,3,5,7,11,13,17,19,23,29,31,37};
unordered_map<int,int>buc;
bool check(int p,int b){
	int d=p-1,tmp=qp(b,d,p);
	if(tmp!=1)return 0;
	while(~d&1){
		if(d>>=1,(tmp=qp(b,d,p))==p-1)return 1;
		else if(tmp!=1)return 0;
	}
	return 1;
}
bool miller(int p){
	if(p<40){
		for(int i=0;i<12;i++)if(pr[i]==p)return 1;
		return 0;
	}
	for(int i=0;i<12;i++)if(!check(p,pr[i]))return 0;
	return 1;
}
int pollard(int x){
	int c=rand()%(x-1)+1,s=0,t=0,d;
	for(int goal=1;;goal<<=1,s=t){
		int prod=1;
		for(int i=1;i<=goal;i++){
			t=((lll)t*t+c)%x;
			prod=(lll)prod*abs(t-s)%x;
			if(i%127==0){if((d=__gcd(prod,x))>1)return d;}
		}
		if((d=__gcd(prod,x))>1)return d;
	}
	return x;
}
void fac(int n,int nk=1){
	if(n<2)return;
	if(miller(n)){buc[n]+=nk;return;}
	int p=0,k=0;
	while((p=pollard(n))==n);
	while(n%p==0)n/=p,k++;
	fac(p,k*nk),fac(n,nk);
}
struct integer {
    int p[4005],len;
    integer(){
        memset(p,0,sizeof(p));
        len=0;
    }
}ans;
integer ti2(integer A){
    integer C;
    for(int i=1,k,f=0;i<=A.len||f;i++){
    	k=f+A.p[i]*2;
    	C.p[i]=k%10,f=k/10;
    	C.len=i;
    }
    return C;
}
void print(integer A){
    if(!A.len){puts("0");return;}
    for(int i=A.len;i>=1;i--)cout<<A.p[i];
    puts("");
}
signed main(){
	srand(20070529);
	for(int T=read();T--;){
		n=read();
		fac(n);
	}
	int ans1=0,ans2=0;
	for(auto t:buc)if(t.second>ans2)ans1=1,ans2=t.second;
	else if(t.second==ans2)ans1++;
	ans.len=1,ans.p[1]=1;
	while(ans1--)ans=ti2(ans);
	cout<<ans2<<'\n';ans.p[1]--;print(ans);
}
posted @ 2022-09-18 10:22  pengyule  阅读(15)  评论(0编辑  收藏  举报