数论一·Miller-Rabin质数测试

#1287 : 数论一·Miller-Rabin质数测试

时间限制: 10000ms

单点时限: 1000ms

内存限制: 256MB

描述

小Hi和小Ho最近突然对密码学产生了兴趣，其中有个叫RSA的公钥密码算法。RSA算法的计算过程中，需要找一些很大的质数。

小Ho：要如何来找出足够大的质数呢？

小Hi：我倒是有一个想法，我们可以先随机一个特别大的初始奇数，然后检查它是不是质数，如果不是就找比它大2的数，一直重复，直到找到一个质数为止。

小Ho：这样好像可行，那我就这么办吧。

过了一会儿，小Ho拿来了一张写满数字的纸条。

小Ho：我用程序随机生成了一些初始数字，但是要求解它们是不是质数太花时间了。

小Hi：你是怎么做的啊？

说着小Hi接过了小Ho的纸条。

小Ho：比如说我要检测数字n是不是质数吧，我就从2开始枚举，一直到sqrt(n)，看能否被n整除。

小Hi：那就对了。你看纸条上很多数字都是在15、16位左右，就算开方之后，也有7、8位的数字。对于这样大一个数字的循环，显然会很花费时间。

小Ho：那有什么更快速的方法么？

小Hi：当然有了，有一种叫做Miller-Rabin质数测试的算法，可以很快的判定一个大数是否是质数。

输入

第1行：1个正整数t，表示数字的个数，10≤t≤50

第2..t+1行：每行1个正整数，第i+1行表示正整数a[i]，2≤a[i]≤10^18

输出

第1..t行：每行1个字符串，若a[i]为质数，第i行输出"Yes"，否则输出"No"

样例输入

3
3
7
9

样例输出

Yes
Yes
No

费马小定理：对于质数p和任意整数a，有a^p
≡ a(mod p)(同余)。反之，若满足a^p
≡ a(mod p)，p也有很大概率为质数。

将两边同时约去一个a，则有a^(p-1)
≡ 1(mod p）

证明：引理1．
　　若a,b,c为任意3个整数,m为正整数,且(m,c)=1,则当ac≡bc(modm)时,有a≡b(modm)
　　证明：ac≡bc(mod
m)可得ac–bc≡0(mod
m)可得(a-b)c≡0(mod
m)因为(m,c)=1即m,c互质,c可以约去,a–
b≡0(mod m)可得a≡b(mod
m)

引理2．
　　设m是一个整数,且m>1,b是一个整数且(m,b)=1.如果a1,a2,a3,a4,…am是模m的一个完全剩余系,则ba[1],ba[2],ba[3],ba[4],…ba[m]也构成模m的一个完全剩余系.
　　证明：若存在2个整数ba和ba[j]同余即ba≡ba[j](mod m),根据引理1则有a≡a[j](mod m).根据完全剩余系的定义可知这是不可能的,因此不存在2个整数ba和ba[j]同余.所以ba[1],ba[2],ba[3],ba[4],…ba[m]构成模m的一个完全剩余系.

构造素数ｐ的完全剩余系   

因为，由引理2可得   

也是p的一个完全剩余系。由完全剩余系的性质，  

即  

易知  ，同余式两边可约去，得到

这样就证明了费马小定理。

二次探测定理：

如果p是奇素数，则x^2≡ 1(mod p)的解为x≡ 1 或x≡ p - 1(mod p)

所以就可以找几个ａ为底，先算一遍费马小定理，然后若p-1为偶数，则用二次探测定理，u=(p-1)/2，检验　x^u
mod p是否为１或　p-1，还要判断ｘ^2modp　要为１，要保证原同余方程是合法的。

还有快速幂的时候可能会乘法溢出，需要搞一个快速加。

 1 LL qmul(LL x, LL y, LL mod)
 2 {
 3   LL ret=0;  
 4   while(y){
 5      if(y&1)
 6      ret=(ret+x)%mod;
 7  　 x=x*2%mod;
 8      y>>=1;
 9   }
10   return ret%mod;  
11 }

可以看成ｙ个ｘ相加，这样可以不打高精。

 1 #include<set>
 2 #include<map>
 3 #include<queue>
 4 #include<stack>
 5 #include<ctime>
 6 #include<cmath>
 7 #include<string>
 8 #include<vector>
 9 #include<cstdio>
10 #include<cstdlib>
11 #include<cstring>
12 #include<iostream>
13 #include<algorithm>
14 #define LL long long
15 using namespace std;
16 int f[13]={2,3,5,7,11,13,17,19,23,29,31,37};
17 LL qmul(LL x, LL y, LL mod)
18 {  
19   LL ret=0;  
20   while(y){  
21     if(y&1)  
22       ret=(ret+x)%mod;  
23     x=x*2%mod;  
24     y>>=1;  
25   }
26   return ret%mod;  
27 }
28 LL POW(LL a,LL p,LL mod){
29   LL ans=1;
30   while(p>0){
31     if(p%2) ans=qmul(ans,a,mod);
32     p>>=1;
33     a=qmul(a,a,mod);
34   }
35   return ans%mod;
36 }
37 bool MR(LL p){
38   for(int i=0;i<5;i++){
39     if(f[i]>=p) break;
40     if(POW(f[i],p-1,p)!=1)return 0;
41     LL pp=p;
42     pp--;
43     while(pp%2==0){
44       pp=pp/2;
45       LL y=POW(f[i],pp,p);
46       if(qmul(y,y,p)==1 && y!=1 && y!=p-1)
47     return 0;
48     }
49   }
50   return 1;
51 }
52 int main()
53 {
54   freopen("!.in","r",stdin);
55   freopen("!.out","w",stdout);
56   int T;scanf("%d",&T);
57   while(T){
58     T--;LL a;
59     scanf("%lld",&a);
60     if(MR(a)) printf("Yes\n");
61     else printf("No\n");
62   }
63   return 0;
64 }

 

转载于:https://www.cnblogs.com/pantakill/p/6629214.html