题面描述:判定一个数P∈[1,2^63-1]∩N是素数么。
按照朴素的判定素数方法,至少也需要O(P^0.5)的,但这道题就是霸气到连这样的时间复杂度都过不了的地步。
实在是不会做了,就学习了传说中的Miller-Rabin素数判定法。
两个引理:
①费马小定理:
设p为质数,且不满足p|a,
则a^(p-1)=a(mod p).
证:
又一个引理,若n与p互质,且a与p互质,则n*a与p互质。
这真的是一个看似很简单的引理,但它却意味着一些看似不那么简单的事情。
设A=(0,p)∩N,则
①对于任意x∈A,x与p互质;
②对于任意x,y∈A,y-x与p互质。
因为假如y-x与p不互质的话,
设y>x,y-x=kp,k∈(0,+无穷)∩N,则y=x+kp,
kp>=p,
∴y∉A,这与题设y∈A相矛盾。
∴y-x与p互质
设B={a*x|x∈A},
则因为a与p互质,
所以对于任意y,x∈A,有Y=a*y∈B,X=a*x∈B,Y-X=a*(y-x)与p互质。
即对于任意x,y∈B,有y-x与p互质。
则C={x(mod p)|x∈B}=A.
因为B中共有p-1项,且均与p互质,所以若其中每项在模p意义下均不相同,则必使得C=A(因为对于任意数x与p互质,模p只可能会产生1~p-1这p-1种可能)。
以下解释为何B中任意两项在模p意义下均不相同:
因为若B中存在两项x=y(mod p),x,y∈B,则p|y-x,但以上我们已证得B中任意两项之差均不能被p整除。
综上,命题得证。
费马小定理是质数的一个性质,但同样存在合数满足此性质。
但我们可以高效地check一个数是否具有此种性质,与快速幂相结合,只需要logP的时间复杂度;但这依然是不够的,因为存在一种叫做伪素数的东西可以很好地卡掉这种算法,而且OI中编数据的人一般对这种伪素数都相当地青睐。
②所以。。我们还需要进一步的检验。
一个与费马小定理完美结合的定理。
模意义下的乘法运算依然遵循四则运算法则,
所以若x^2=1,则x=±1,在模运算中,就是:
若x^2=1(mod p),则x=1或p-1.
说这玩意儿高明,究竟高明在哪儿呢?
①它与费马小定理完美结合,完全相同的运算方式。
②在最好情况下,它可以实现log^2级的高效check。
所以呢,我们的做法就是随机选几个数,然后判断它是否符合费马小定理,若符合再将指数除以2不断进行二次检验;直到检验出其不是质数或结果等于-1为止。
综上所述,我们很容易就可以发现这本质上其实是一种hash的思想,通过两个异常巧妙的结合。
这就是Miller-Rabin素数判定法的全部奥妙所在了;据说每一次检验会有1/4的可能出错,所以如果是对于2^63-1这样的数据范围的话,选33个数以上应该是比较合理的。。
时间复杂度为,在下面的代码中选取了10个数,并把高精度乘法改成了写起来更为简单的二分求积,所以
#include<iostream>
using namespace std;
#include<cstdlib>
#include<cstdio>
typedef unsigned long long LLU;
LLU N;
inline LLU multi(LLU x,LLU y){
LLU a=x,b=y,t=x,ans=0;
for(;b;b>>=1){
if(b&1)ans=(ans+t)%N;
t=(t<<1)%N;
}
//cout<<x<<"*"<<y<<"="<<ans<<endl;
return ans;
}
inline bool Prime(LLU A,LLU tmp){
//cout<<"Prime:"<<A<<" "<<tmp<<endl;
LLU a=A,b=tmp,t=a,ans=1;
for(;b;b>>=1){
if(b&1)ans=multi(ans,t);
t=multi(t,t);
}
//cout<<A<<"^"<<tmp<<"="<<ans<<endl;
if(ans==1)return 0;
if(ans==N-1){
if(tmp==N-1){
printf("No");
exit(0);
}
return 1;
}
else{
printf("No");
exit(0);
}
}
int main(){
LLU tmp,a,b,t,ans,j;
int i,prime[]={10,35,77,535,71497,2,3,5,7,11,3161};
cin>>N;
if(N==1){
printf("No");
return 0;
}
for(i=prime[0];i--;)
if(N%prime[i]&&prime[i]%N){
j=N-1<<1;
do{
j>>=1;
if(Prime(prime[i],j))
break;
}while(~j&1);
}
printf("Yes");
}