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);
}