素数判定Miller_Rabin 算法详解
上次说好的要把素数判定和大数分解（见另一篇博文）的快速随机化算法解决了，于是乎今天就来解决，不得不说理解起来真的有困难。我只能大概的将思路理一下，若有错漏还请担待。
首先相信有一些数学和编程经验的读者应该知道，最简单直观简单的素数判定方法就是试除法。对于判断数n是否是素数，我们从2开始一直到sqrt(n)。如果找到一个因子则判断n不是素数，否则是素数。代码如下：
- bool isPrime( long long n )
- {
- for(long long i = 2; i*i <= n; i++)
- {
- if(n%i == 0) return false;
- }
- return true;
- }
bool isPrime( long long n )
{
for(long long i = 2; i*i <= n; i++)
{
if(n%i == 0) return false;
}
return true;
}
所以随着学习的深入，我们了解到了素数筛法，即从2开始，2的倍数肯定不是素数，再向右扫描，如果扫描到素数，则重复之前的过程，剔除之后的部分合数（准确的说是关于当前质数的倍数），如果扫描到合数则跳过（表示前面已经更新过这个数不是素数）。然后都扫描一遍即可把1~n的素数求解出来。这个算法的复杂度略高于O(n)。素数筛代码如下：
- #include<iostream>
- #include<cstring>
- using namespace std;
- const int MAXN = 500000;
- bool isprime[MAXN];
- int prime[MAXN];
- int cnt = 0;//保存素数个数
- void getPrime()
- {
- for(int i = 1; i < MAXN; i++)
- isprime[i] = true; //先假设所有数是素数，后面逐个扫描更新
- for(int i = 2; i < MAXN; i++) //扫一遍
- {
- if(!isprime[i]) continue; //如果不是素数，则不往后面更新
- prime[cnt++] = i;
- for(int j = 2 * i; j < MAXN; j += i)
- isprime[j] = false;
- }
- }
- int main()
- {
- getPrime();
- for(int i = 0; i < cnt; i++)
- cout << prime[i] << endl;
- }
#include<iostream>
#include<cstring>
using namespace std;
const int MAXN = 500000;
bool isprime[MAXN];
int prime[MAXN];
int cnt = 0;//保存素数个数
void getPrime()
{
for(int i = 1; i < MAXN; i++)
isprime[i] = true; //先假设所有数是素数，后面逐个扫描更新
for(int i = 2; i < MAXN; i++) //扫一遍
{
if(!isprime[i]) continue; //如果不是素数，则不往后面更新
prime[cnt++] = i;
for(int j = 2 * i; j < MAXN; j += i)
isprime[j] = false;
}
}
int main()
{
getPrime();
for(int i = 0; i < cnt; i++)
cout << prime[i] << endl;
}
但是这个算法的弊端在于，为了判断一个大数是否是素数必须从从头开始扫描，而且空间代价也受不了，故不能离散的判断。
---------------------------------------------------------------------------------------------------------------------------------------------
于是进入我们今天的主题，Miller_rabin算法，优势可以单独判断一个大数是否素数。缺点他是一个不保证正确的算法，我们只能通过多次执行算法让这个错误的概率很小，不过幸运的是通常来看它的错误概率可以小到忽略不计。
Miller_rabin算法描述
算法的理论基础：
1. Fermat定理：若n是奇素数，a是任意正整数(1≤ a≤ n−1)，则 a^(n-1) ≡ 1 mod n。
2. 推演自Fermat定理（具体过程我没看懂，Orz）, 如果n是一个奇素数，将n−1表示成2^s*r的形式，r是奇数，a与n是互素的任何随机整数，那么a^r ≡ 1 mod n或者对某个j (0 ≤ j≤ s−1, j∈Z) 等式a^(2jr) ≡ −1 mod n 成立。
聪明的读者已经发现上述的定理是一个数是素数的必要条件而非充分条件，所以这就是这个算法的不精确的原有，对于一些特定的检验算子a，存在一些合数也满足上述定理。所以我们要多找几个a反复检验，这样能让错误率大大降低。
那么我们按照上述的定理2，首先重复n次实验。对于每一次实验，随机取检验算子a，带入定理2进行检验，看看在算子a下，n能否满足
a^r ≡ 1 mod n或者对某个j (0 ≤ j≤ s−1, j∈Z) 等式a^(2jr) ≡ −1 mod n **
如果任意一次实验不满足，则判定不是素数，如果都满足，可近似可以认为是素数（错误率极小）。
代码实现如下：（代码中的q_mul和q_pow表示快速乘法和快速幂运算，具体讲解请参考另一篇博文）
- #include <iostream>
- #include <cstdio>
- #include <algorithm>
- #include <cmath>
- #include <cstring>
- #include <map>
- using namespace std;
- const int times = 20;
- int number = 0;
- map<long long, int>m;
- long long Random( long long n ) //生成[ 0 , n ]的随机数
- {
- return ((double)rand( ) / RAND_MAX*n + 0.5);
- }
- long long q_mul( long long a, long long b, long long mod ) //快速计算 (a*b) % mod
- {
- long long ans = 0;
- while(b)
- {
- if(b & 1)
- {
- b--;
- ans =(ans+ a)%mod;
- }
- b /= 2;
- a = (a + a) % mod;
- }
- return ans;
- }
- long long q_pow( long long a, long long b, long long mod ) //快速计算 (a^b) % mod
- {
- long long ans = 1;
- while(b)
- {
- if(b & 1)
- {
- ans = q_mul( ans, a, mod );
- }
- b /= 2;
- a = q_mul( a, a, mod );
- }
- return ans;
- }
- bool witness( long long a, long long n )//miller_rabin算法的精华
- {//用检验算子a来检验n是不是素数
- long long tem = n - 1;
- int j = 0;
- while(tem % 2 == 0)
- {
- tem /= 2;
- j++;
- }
- //将n-1拆分为a^r * s
- long long x = q_pow( a, tem, n ); //得到a^r mod n
- if(x == 1 || x == n - 1) return true; //余数为1则为素数
- while(j--) //否则试验条件2看是否有满足的 j
- {
- x = q_mul( x, x, n );
- if(x == n - 1) return true;
- }
- return false;
- }
- bool miller_rabin( long long n ) //检验n是否是素数
- {
- if(n == 2)
- return true;
- if(n < 2 || n % 2 == 0)
- return false; //如果是2则是素数，如果<2或者是>2的偶数则不是素数
- for(int i = 1; i <= times; i++) //做times次随机检验
- {
- long long a = Random( n - 2 ) + 1; //得到随机检验算子 a
- if(!witness( a, n )) //用a检验n是否是素数
- return false;
- }
- return true;
- }
- int main( )
- {
- long long tar;
- while(cin >> tar)
- {
- if(miller_rabin( tar )) //检验tar是不是素数
- cout << "Yes, Prime!" << endl;
- else
- cout << "No, not prime.." << endl;
- }
- return 0;
- }
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <map>
using namespace std;
const int times = 20;
int number = 0;
map<long long, int>m;
long long Random( long long n ) //生成[ 0 , n ]的随机数
{
return ((double)rand( ) / RAND_MAX*n + 0.5);
}
long long q_mul( long long a, long long b, long long mod ) //快速计算 (a*b) % mod
{
long long ans = 0;
while(b)
{
if(b & 1)
{
b--;
ans =(ans+ a)%mod;
}
b /= 2;
a = (a + a) % mod;
}
return ans;
}
long long q_pow( long long a, long long b, long long mod ) //快速计算 (a^b) % mod
{
long long ans = 1;
while(b)
{
if(b & 1)
{
ans = q_mul( ans, a, mod );
}
b /= 2;
a = q_mul( a, a, mod );
}
return ans;
}
bool witness( long long a, long long n )//miller_rabin算法的精华
{//用检验算子a来检验n是不是素数
long long tem = n - 1;
int j = 0;
while(tem % 2 == 0)
{
tem /= 2;
j++;
}
//将n-1拆分为a^r * s
long long x = q_pow( a, tem, n ); //得到a^r mod n
if(x == 1 || x == n - 1) return true; //余数为1则为素数
while(j--) //否则试验条件2看是否有满足的 j
{
x = q_mul( x, x, n );
if(x == n - 1) return true;
}
return false;
}
bool miller_rabin( long long n ) //检验n是否是素数
{
if(n == 2)
return true;
if(n < 2 || n % 2 == 0)
return false; //如果是2则是素数，如果<2或者是>2的偶数则不是素数
for(int i = 1; i <= times; i++) //做times次随机检验
{
long long a = Random( n - 2 ) + 1; //得到随机检验算子 a
if(!witness( a, n )) //用a检验n是否是素数
return false;
}
return true;
}
int main( )
{
long long tar;
while(cin >> tar)
{
if(miller_rabin( tar )) //检验tar是不是素数
cout << "Yes, Prime!" << endl;
else
cout << "No, not prime.." << endl;
}
return 0;
}
5886
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
