费马小定理(即若P为质数,则\(A^P\equiv A \pmod{P}\))。
欧几里得算法(GCD)。
快速幂,龟速乘。
素性测试是OI中一个十分重要的事,在数学毒瘤题中有着举足轻重的地位。
常见的素性测试如下:
int check(int N){
for(int i=2;i*i<=N;i++)
if(N%i==0)return 0;
return 1;
}
以上是一个\(O(\sqrt{N})\)的算法,虽然不优,但在绝大多数情况下是可以的。
但是,假如\(N\)的范围达到了\(1e18\),以上算法很明显就不行了,我们得考虑更优的算法。
引入Miller-Rabin算法。
定理1:如果\(P\)为一个大于2的素数,那么方程\(X^2\equiv1 \pmod{P}\)的解只有1或者-1。
(是真得水的) 证明如下:
由\(X^2\equiv1 \pmod{P}\),得\((X^2-1)\equiv0 \pmod{P}\),
即\({(X+1)(X-1)}\equiv0 \pmod{P}\)。
则\(P \mid (X+1)\)或\(P\mid(X-1)\),
即\((X+1)\equiv0 \pmod{P}\)或\((X-1)\equiv0 \pmod{P}\)。
则\(X\equiv-1 \pmod{P}\)或\(X\equiv1 \pmod{P}\),即证。
(注:由于P>2,不能同时满足\(X\equiv-1 \pmod{P}\)与\(X\equiv1 \pmod{P}\)。)
定理1的逆否定理即为:
如果有\(P>2\)且\(X^2\equiv1 \pmod{P}\)且\(X\not\in\{1,-1\}\),则\(P\)不为质数。
证明如下:
假设\(P\)为质数,则\(P\)满足定理1,即方程\(X^2\equiv1 \pmod{P}\)的解只有1或者-1。
与条件矛盾,故\(P\)不为质数。
基于以上定理,我们可以开始了。
假如我们要验证\(P(P>2)\)是否是一个质数。
先设\(P\)是一个素数。
那么对于所有满足\(A^2\equiv1 \pmod{P}\)的整数\(A\),都有\(A\in\{1,-1\}\)。
但如果我们枚举\(A\)来暴力判断,则复杂度又上升到了\(O(N)\)。
(不妨令\(A<P\),以下\(A\)均小于\(P\))
考虑做判断时,如何快速的找到满足\(A^2\equiv1 \pmod{P}\)的整数\(A\)。
我们考虑将\((P-1)\)分解成\(2^S·D\)的形式,其中\(D\)为奇数。
(由于\(P>2\)且\(P\)为素数,则\((P-1)\)为偶数,即\(S\ge1\))
则根据费马小定理,有
\[A^{P-1}\equiv 1\pmod{P}
\]
\[A^{2^S·D}\equiv 1\pmod{P}
\]
\[({A^{2^{S-1}·D}})^2\equiv 1\pmod{P}
\]
则此时\(A^{2^{S-1}·D}\)又变成了在\(A^2\equiv1 \pmod{P}\)的新的\(A\)值。
不妨设\(A\)值数组为\(\{A_0,A_2,A_3,A_4….A_{S-1}\}\),其中\(A_i=A^{2^i·D}\)。
(注:不难发现对于任意满足\(0<i<S\)的\(i\),都有\(A_i=2·A_{i-1}\))
这样的话对于每个初始值\(A\),我们可以知道所构成的\(A\)数组值。
此时我们可以考虑从\(A_0\)走到\(A_{S-1}\)(即从\(A^D\)到\(A^{2^{S-1}·D}\))。
①考虑初始情况的\(A_0\)(即为\(A^D\)),
它如果满足\(A_0\equiv1 \pmod{P}\)或\(A_0\equiv-1 \pmod{P}\)。
则对于任意\(A_i\)(\(0<i< S\)),都有\(A_i\equiv1 \pmod{P}\),即通过本次测试。
(注:通过测试只能说明它可能为素数,而没通过则说明它一定不是素数。)
②考虑过程中的情况:设有一满足\(0<t<S\)的\(t\)。
若有\(A_t\equiv-1\pmod{P}\),则说明
对于任意满足\(t<i<S\)的\(i\),都有\(A_i\equiv1\pmod{P}\),即通过此次测试。
若有\(A_t\equiv1\pmod{P}\),则说明\({A_{t-1}}^2\equiv1\pmod{P}\),
而\(A_{t-1}\)又不满足\(A_{t-1}\equiv-1\pmod{P}\)或\(A_{t-1}\equiv-1\pmod{P}\),
因为如果\(A_{t-1}\)满足以上条件,则在上一次或\(A_0\)时就已经判了。
故存在一个\(A\),使得\(A^2\equiv1\pmod{P}\)且\(A\not\in\{1,-1\}\)(定理1的延伸)。
故未通过此次测试。
注:若\(P(P>2)\)为素数: 则\(A^2\equiv1\pmod{P}\)只有\(A\)满足
\(A\equiv1\pmod{P}\)或\(A\equiv-1\pmod{P}\)时满足(由定理1得)。
但\(A^2\equiv-1\pmod{P}\)对于\(A\)的取值却没有特殊性。
③考虑结尾情况:
若对于这一组\(A\)都没有一个\(A\)使得\(A^2\equiv1\pmod{P}\),
则在判\(A_{S-1}\)时,实际上判的是\({({A^{2^{S-1}·D}})}^2\not\equiv1\pmod{P}\)
即\(A^{P-1}\not\equiv1\pmod{P}\),故\(P\)不满足费马小定理,即\(P\)未通过测试。
为保证算法的正确性与高效性,
我们考虑随机化许多个整数\(A\)来多次判断使得\(P\)为素数的概率尽量大。
那么我们要随机多少次才能使得\(P\)素数的概率能让我们接受呢?
根据大佬们的研究表明,Prime是一个RP问题(通过大量数据来算概率),
而对于每次测试的期望失败概率是比\(\frac{1}{4}\)要小的。(通过试验得到的结论)
即我们做K次成功的测试然后判错\(P\)的概率即为\(4^{-K}\)。
这个失败率在代码中也是可以接受的。
(注:经过试验,当A值取\(\{2,3,5,7,11,13,17,19,23,29\}\)时,可以过long long
级别的询问)
int MillerTest(LL N){//素数测试
if(N<=1||N==4)return 0;//特判N比较小情况
if(N<=3)return 1;
LL D=N-1;
while(D%2==0)D/=2;
for(int i=1;i<=PK;i++){
if(Prime[i]>=N)return 1;
int Ok=0;
LL A=Prime[i];
LL x=quick_Pow(A,D,N);
if(x==1||x==N-1)
continue;
while(D!=N-1){
x=quick_Mul(x,x,N);//龟速乘防爆long long
D*=2;
if(x==1)return 0;
if(x==N-1){Ok=1;break;}
}
if(!Ok)return 0;
}
return 1;
}
质因数分解是OI中一个十分重要的事,在数学毒瘤题中有着举足轻重的地位。
常见的质因数分解如下:
void divide(int N){
for(int i=2;i*i<=N;i++){
int cnt=0;
while(N%i==0)cnt++,N/=i;
if(cnt)printf("%d %d\n",i,cnt);
}
if(N!=1)printf("%d 1\n",N);
}
以上是一个\(O(\sqrt{N})\)的算法,虽然不优,但在绝大多数情况下是可以的。
但是,假如\(N\)的范围达到了\(1e18\),以上算法很明显就不行了,我们得考虑更优的算法。
引入Pollard-Rho算法。
说起Pollard-Rho算法,我们可能不太熟悉。
但提到随机化,大家应该多少都有些了解。
//随机化找N的一个因数
int i=0;do{
i=rand()%(N-2)+2;
}while(N%i);
printf("I found %d \n",i);
我们不难发现:随机化一次找到某个因数的概率十分小。
当\(N=P·Q\)时(\(P,Q\)均为素数),成功的概率为\(\frac{2}{N-1}\)。
期望次数为\(O(N)\)级别,(甚至不如暴力)。
但是,如果我们加上一些优化呢?
考虑每次随机化。
其实对于随机的某个数\(i\),只要\(gcd(i,N)\not=1\)就可以找到一个因数。
(时间复杂度仿佛还是\(O(N)\)….)
我们考虑用生日悖论来优化随机化。
尽管听起来特别玄学,但是此做法还是有一定科学依据的。
生日悖论,指如果一个房间里有23个或23个以上的人,那么至少有两个人的生日相同的概率要大于50%。这就意味着在一个典型的标准小学班级(30人)中,存在两人生日相同的可能性更高。对于60或者更多的人,这种概率要大于99%。
以上是生日悖论的一个例子,其函数图像大致如下。
考虑换一种思考方式,我们随机化\(N\)个数出来,得到\(N^2\)个差值。
然后我们用这些差值来执行优化1,时间复杂度将明显降低。
考虑生日相同的情况,其可以被描述为某两个人的生日值差值为0。
从以上描述差值为0的图像中,可以明显发现选取的人数\(N\)越大,则概率越大。
则可以发现\(N\)的升大对于所有\(K\)的概率均有提升。
(易发现\(N^2\)对于每个\(K\)的期望值为总值域\(W\)的根号\(O(\sqrt{W})\))
然后,我们使\(N\)值固定。
我们设\(P(K)\)表示以\(K\)为差值的概率。
则感性发现:\(P(1)>P(0),P(W-1)>P(W)\)
感性推出:实际的\(P\)图像其实是在前方激增,中间接近1,后方再次猛降。
故我们在判断时,得特判差值出现在前面的情况。
(注:\(N\)越大时,\(P\)图像的前方激增斜率便越大。)
应用优化1与优化2,考虑如何找到\(N^2\)个差值,
(可以认为生日悖论中的某个差值\(K\)是我们所需要的答案)
我们考虑随机出两个数\(X\),\(Y\),然后用\(abs(X-Y)\)作为差值。
然后令\(f(X)=(X^2+P) \mod W,X=f(X),Y=f(Y)\)。
这样的话,我们只用随机化两个初值,一个常数\(P\),就可以了。
但是易发现这样可能重复,即走进一个无解的循环里。
所以我们让\(X,Y\)的初值相等,每次随机的时候令\(Y=f(f(Y))\),
即出发点一样,\(Y\)的速度为\(X\)的两倍,若有无解循环,
则总有一次\(X,Y\)会相遇(复杂度至多为经过的点数量),此时我们就可以直接退出了。
(注:\(f(X)=(X^2+P) \mod W\)是一个先人传下的神奇函数)
(其神奇之处在于玄学般地大幅提升了代码效率)
其实随机这些差值的本质还是随机化,依旧是\(N^2\)个差值,但复杂度却玄学降低了。
???
故这个算法的复杂度十分看RP,
在一些数据的复杂度甚至超过了\(O(\sqrt{W}\)),
不过能解决一些long long
的数据还是很强大了。
结合Miller-Rabin与Pollard-Rho算法。
我们对于一个整数\(N\),先用Miller-Rabin判断其是否为质数。
若是,则上传N。
若不是,则用Pollard-Rho找一个因数出来,
将其分解成更小的两个数继续递推。
//该代码只是找到了N的一个因数,可能并不是素数
LL PollardRho(LL N){//找因数
if(N==1)return -1;
for(int i=1;i<=PK;i++)
if(N%Prime[i]==0)
return Prime[i];//特判出奇迹
LL x=(rand()%(N-2))+2,y=x;
LL P=(rand()%(N-1))+1,D=1;
while(D==1){
x=(quick_Mul(x,x,N)+P)%N;
y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
D=gcd(abs(x-y),N);//x=y时,D值为N
}
return D;
}
板题:Prime Test
大意:找到\(N\)的最小的质因数。
#include<ctime>
#include<cstdio>
#include<algorithm>
using namespace std;
#define LL long long
const int MAXK=15;
const long long INF=1e18;
int PK=10,Prime[MAXK]={0,2,3,5,7,11,13,17,19,23,29};
//此行数对于long long范围内的素性测试已够用了
long long N;
LL Add(LL x,LL y,LL p){
if(x+y<p)return x+y;
return x+y-p;
}
LL gcd(LL x,LL y){//最大公因数
if(x==0)return y;
return gcd(y%x,x);
}
LL quick_Mul(LL x,LL y,LL p){//龟速乘
x%=p;y%=p;
LL ret=0;
while(y){
if(y&1)ret=Add(ret,x,p);
x=Add(x,x,p);
y>>=1;
}
return ret;
}
LL quick_Pow(LL x,LL y,LL p){//快速幂
LL ret=1;
while(y){
if(y&1)ret=quick_Mul(ret,x,p);
x=quick_Mul(x,x,p);
y>>=1;
}
return ret;
}
int MillerTest(LL N){//素数测试
if(N<=1||N==4)return 0;
if(N<=3)return 1;
LL D=N-1;
while(D%2==0)D/=2;
for(int i=1;i<=PK;i++){
if(Prime[i]>=N)return 1;
int Ok=0;
LL A=Prime[i];
LL x=quick_Pow(A,D,N);
if(x==1||x==N-1)
continue;
while(D!=N-1){
x=quick_Mul(x,x,N);
D*=2;
if(x==1)return 0;
if(x==N-1){Ok=1;break;}
}
if(!Ok)return 0;
}
return 1;
}
LL PollardRho(LL N){//找因数
if(N==1)return -1;
for(int i=1;i<=PK;i++)
if(N%Prime[i]==0)
return Prime[i];//特判出奇迹
LL x=(rand()%(N-2))+2,y=x;
LL P=(rand()%(N-1))+1,D=1;
while(D==1){
x=(quick_Mul(x,x,N)+P)%N;
y=(quick_Mul(y,y,N)+P)%N,y=(quick_Mul(y,y,N)+P)%N;
D=gcd(abs(x-y),N);
}
return D;
}
LL Min_Div(LL N){//最小质因数
if(N==1)return INF;
if(MillerTest(N))return N;
LL X=PollardRho(N);
return min(Min_Div(X),Min_Div(N/X));
}
int main(){
//srand(time(NULL));
while(~scanf("%lld",&N)){
printf("%lld\n",Min_Div(N));
}
}
手机扫一扫
移动阅读更方便
你可能感兴趣的文章