从龟速乘到 $Miller-Rabin$ 算法(数论算法总结)
发现自己竟然菜到不太会龟速乘,所以把 \(Miller-Rabin\) 算法所需要用到的算法全学了一遍……
龟速乘
龟速乘是一种 \(O(\log n)\) 的乘法计算方法。
考虑有时普通乘法取模会爆 \(long\ long\),因此我们考虑用类似快速幂的方式进行乘法运算。
int mul(int x,int y,int c){
x%=c,y%=c;
int re=0;
while(y){
if(y&1) re+=x,re-=(re>=c)?c:0;
x+=x,x-=(x>=c)?c:0,y>>=1;
}return re;
}
快速乘
快速乘是一种 \(O(1)\) 的乘法计算方法。
发现 \(long\ double\) 可以存储 \(10^{300}\) 的数,所以考虑使用 \(long\ double\) 进行乘法运算。我还没有遇到过出锅的情况。
int mul(int a,int b,int p){
return (a*b-(int)((long double)a/p*b)*p+p)%p;
}
\(Miller-Rabin\) 算法
哈,这玩意错的跟对的一样。
和快速幂求逆元一样,我们考虑费马小定理。
我们可以得到:
设 \(f(x)=a^{x-1}\bmod x\),则当 \(f(p)\ne 1\) 时,\(p\) 一定是合数;当 \(p\) 是质数时,\(f(p)=1\)。
但是正确率不是人类所能接受的,所以考虑优化。
考虑下面这个性质:
若 \(b^2\bmod p=1\),\(p\) 为质数,则 \(b-1|p\) 或 \(b+1|p\)。
这样我们就可以进行二次探测,正确率大大提高。
我们设用于测试的集合为 \(A\),则 \(A=\{2,3,5,7,11,13,17,19,23\}\) 时,对于 \(x\le 10^{18}\) 的情况,都可以正确判断 \(x\) 的素性。
int mr[9]={2,3,5,7,11,13,17,19,23};
int mul(int a,int b,int p){
return (a*b-(int)((long double)a/p*b)*p+p)%p;
}int qpow(int x,int y,int c){
int re=1;
while(y){
if(y&1) re=mul(re,x,c);
x=mul(x,x,c),y>>=1;
}return re;
}int check(int x,int md){
int c=md-1,mid=qpow(x,c,md);
if(mid!=1) return 0;
while(!(c&1)&&mid==1)
c>>=1,mid=qpow(x,c,md);
return (mid==1||mid==md-1);
}int miller_rabin(int x){
if(x<2) return 0;
if(x<=23){
for(int i=0;i<9;i++)
if(mr[i]==x) return 1;
return 0;
}for(int i=0;i<9;i++)
if(!check(mr[i],x)) return 0;
return 1;
}
\(ps.Pollard\ rho\) 不太懂,贴份模板题代码,先咕咕了。
#include<bits/stdc++.h>
#define int long long
using namespace std;
int k,n;
int mr[9]={2,3,5,7,11,13,17,19,23};
int mul(int a,int b,int p){
return (a*b-(int)((long double)a/p*b)*p+p)%p;
}int qpow(int x,int y,int c){
int re=1;
while(y){
if(y&1) re=mul(re,x,c);
x=mul(x,x,c),y>>=1;
}return re;
}int check(int x,int md){
int c=md-1,mid=qpow(x,c,md);
if(mid!=1) return 0;
while(!(c&1)&&mid==1)
c>>=1,mid=qpow(x,c,md);
return (mid==1||mid==md-1);
}int miller_rabin(int x){
if(x<2) return 0;
if(x<=23){
for(int i=0;i<9;i++)
if(mr[i]==x) return 1;
return 0;
}for(int i=0;i<9;i++)
if(!check(mr[i],x)) return 0;
return 1;
}int gcd(int x,int y){
return (!y)?x:gcd(y,x%y);
}int pr(int x){
int s=0,t=0,val=1;
int c=rand()%(x-1)+1;
for(int gl=1;;gl*=2,s=t,val=1){
for(int st=1;st<=gl;st++){
t=(mul(t,t,x)+c)%x;
val=mul(val,abs(t-s),x);
if(st%127==0){
int d=gcd(val,x);
if(d>1) return d;
}
}int d=gcd(val,x);
if(d>1) return d;
}
}void pol_rho(int x,int &mx){
if(x<=mx||x<2) return;
if(miller_rabin(x))
return mx=max(mx,x),void();
int p=x;while(p>=x) p=pr(x);
while(x%p==0) x/=p;
pol_rho(x,mx),pol_rho(p,mx);
}signed main(){
ios::sync_with_stdio(0);
cin.tie(0),cout.tie(0);
srand(time(0));
cin>>k;
while(k--){
cin>>n;int ans=0;
pol_rho(n,ans);
if(ans==n) cout<<"Prime\n";
else cout<<ans<<"\n";
}return 0;
}