#1287 : 数论一·Miller-Rabin质数测试
描述
小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^p ≡ a(mod p)(同余)。反之，若满足a^p ≡ a(mod p)，p也有很大概率为质数。 将两边同时约去一个a，则有a^(p-1) ≡ 1(mod p);
因此，可以用这个定理来测试。任选一个小于p的正整数，检查a^(p-1) mod p是否等于1，若不等于，则100%不是素数。若只测一次，有1/4的错误率，因此，多拿几个a测试，便可大大提高正确性，该测试被称为Fermat测试。
但是：Fermat测试不一定是准确的，有可能出现把合数误判为质数的情况。Miller和Rabin在Fermat测试上，建立了Miller-Rabin质数测试算法。
与Fermat测试相比，增加了一个二次探测定理：如果p是奇素数，则 x^2 ≡ 1(mod p)的解为 x ≡ 1 或 x ≡ p – 1(mod p)。
即，若一个小于p的正整数x，满足x^2≡1mod(p) ,则解一定是x==1或x==p-1,若满足先前的条件，但解却不是这个，则说明这个肯定不是素数，（至于为什么，我也不知道）。举个栗子。假设n=341，我们选取的a=2。则第一次测试时，2^340 mod 341=1。由于340是偶数，因此我们检查2^170，得到2^170 mod 341=1，满足二次探测定理。同时由于170还是偶数，因此我们进一步检查2^85 mod 341=32。此时不满足二次探测定理，因此可以判定341不为质数。
于是，遍得到了一下算法：
for(int i=0;i<5;i++) { ll np=n; ll a=ss[i]; if(a>=n) break; if(calmod(a,n-1,n)!=1) return 0; np--; while(np%2==0) { np>>=1; ll y=calmod(a,np,n); if(mutiply(y,y,n)==1&&y!=1&&y!=n-1) return 0;// 同时满足,x^2 mod p ==1,且y==1||y==n-1,若不满足后面，则克判断非素数！！ } }
View Code
另外，快速幂里面相乘需要用到快速加（原理与快速幂一样），防止抱longlong；
1 #define ll long long 2 #include<algorithm> 3 #include<iostream> 4 #include<iomanip> 5 #include<cstring> 6 #include<cstdlib> 7 #include<cstdio> 8 #include<queue> 9 #include<ctime> 10 #include<cmath> 11 #include<stack> 12 #include<map> 13 #include<set> 14 using namespace std; 15 ll ss[5]={2,3,5,7,11}; 16 ll mutiply(ll a,ll b,ll mod) { 17 ll res=0; 18 while(b){ 19 if(b&1) res=(res+a)%mod; 20 b>>=1; 21 a=a*2%mod; 22 } 23 return res%mod; 24 } 25 ll calmod(ll a,ll n,ll mod) { 26 ll res=1; 27 while(n) { 28 if(n&1) res=mutiply(res,a,mod); 29 n>>=1; 30 a=mutiply(a,a,mod); 31 } 32 return res%mod; 33 } 34 bool MRtest(ll n) { 35 if(n==2) return 1; 36 for(int i=0;i<5;i++) { 37 ll np=n; 38 ll a=ss[i]; 39 if(a>=n) break; 40 if(calmod(a,n-1,n)!=1) return 0; 41 np--; 42 while(np%2==0) { 43 np>>=1; 44 ll y=calmod(a,np,n); 45 if(mutiply(y,y,n)==1&&y!=1&&y!=n-1) return 0;// 同时满足,x^2 mod p ==1,且y==1||y==n-1,若不满足后面，则克判断非素数！！ 46 } 47 } 48 return 1; 49 } 50 int main() {// please remember :infer other things from one fact 51 freopen("Miller-Rabin.in","r",stdin); 52 freopen("Miller-Rabin.out","w",stdout); 53 int t;scanf("%d",&t); 54 while(t--) { 55 ll n; 56 scanf("%lld",&n); 57 if(MRtest(n)) puts("Yes"); 58 else puts("No"); 59 } 60 return 0; 61 }
View Code