判断质数 _ Miller Rabin算法
算法作用
判断一个数是否是素数
算法依据
费马小定理
如果P是素数，且整数a不是p的倍数有：
\[a^{p-1}\equiv 1 \pmod P
\]
费马定理只是n是素数的必要条件。即费马定理不成立，n一定是合数；费马定理成立，n可能是素数。
点击查看证明
性质 1：$p-1$个整数 $a,2a,3a,...(p-1)a$中没有一个是 $p$的倍数
性质 2：$a,2a,3a,...(p-1)a$中没有任何两个同余与模 $p$的
所以 $a,2a,3a,...(p-1)a$对模 $p$的同余既不为零，也没有两个同余相同
因此，这 $p-1$个数模 $p$的同余一定是 $a,2a,3a,...(p-1)a$的某一种排列
即 $a*2a*3a*...*(p-1)a \equiv {1*2*3*...*(p-1)} \pmod p$
化简为
$a^{p-1}*(p-1)! \equiv {p-1}! \pmod p$
根据威尔逊定理可知 $(p-1)!$与 $p$互质，所以同时约去 $(p-1)!$
即得到 $a^{p-1}\equiv 1 \pmod P$
二次探测定理
若 \(p\)为素数，对于：
\[a^{2}\equiv 1 \pmod P
\]
求解的解只有两种可能：a = 1 或者 a = p-1
点击查看证明
$a^{2}\equiv 1 \pmod P$
$a^{2}-1\equiv 0 \pmod P$
$(a+1)*(a-1)\equiv 0 \pmod P$
那么
$(a+1)\equiv 0 \pmod P$
或者
$(a-1)\equiv 0 \pmod P$
(此处可根据唯一分解定理证明)
即
$a\equiv \pm 1 \pmod P$
$a^{2}\equiv 1 \pmod P$
$a^{2}-1\equiv 0 \pmod P$
$(a+1)*(a-1)\equiv 0 \pmod P$
那么
$(a+1)\equiv 0 \pmod P$
或者
$(a-1)\equiv 0 \pmod P$
(此处可根据唯一分解定理证明)
即
$a\equiv \pm 1 \pmod P$
推理定理
Miller Rabin 算法的过程
假设需要判断的数是 \(p\)
我们把 \(p-1\)分解为 \(2^k*t\)的形式
当 \(p\)是素数，有 \(a ^ {2^k * t} \equiv 1 \pmod p\)
然后随机选择一个数 \(a\)，计算出 \(a^t \pmod p\)
让其不断的乘 \(a\)，同时结合二次探测定理进行判断
如果
\[自乘后的数\pmod p = 1
\]
那么这个数就是合数 (违背了二次探测定理)
假如p是素数，那么自乘后的数应该== 1或n-1
这样乘 \(k\)次，最后得到的数就是 \(a^{p-1}\)
那么如果最后计算出的数不为 \(1\)，这个数也是合数 (费马小定理)
正确性判断
老祖宗告诉我们，若 \(p\)通过一次测试，则 \(p\)不是素数的概率为 \(25\)%
那么经过 \(t\)轮测试,\(p\)不是素数的概率为 \(\dfrac {1}{4^{t}}\)
我习惯用 \(2,3,5,7,11,13,17,19\)这几个数进行判断
在信息学范围内出错率为 \(0\)%(不带高精)
算法流程
- 首先对于要判断的数 n，先判断是不是 2，是就直接返回 true。
- 判断是不是小于 2 或者为偶数，是就返回 false。
- 令 n-1 = m*2t，求出 u 和 t。
- 利用 rand 随机找一个 a， a 属于 [1， n-1]
- 令 x = (am)%n
- 然后进行 t 次循环，为什么要循环 t 次呢，上面证明是从 2t 倒着推到 20，我们写程序完全可以正着来，上面的 x 其实就是
- 那么在循环中每次都进行 y =（x*x）%n， x = y 的操作，最终会变为 an-1，也就是
又因为 x 在循环时是一定小于 n 的，所以可以用二次探测定理
运算过程中如果 (x2)%n = 1, 也就是 y = 1，假如 n 是素数，那么 x == 1 或者 x == n-1，否则 n 就一定不是素数，直接返回 false。 - 整个循环结束后，程序运行到最后 x = an-1, 根据费马小定理，如果 x 还是不等于 1，那么肯定不是素数，直接返回 false。
- Miller-Rabin 算法正确率为 75%，所以要进行大约 4~8 次来提高正确率。
代码
#include <iostream>
#include <cstdlib>
using namespace std;
typedef long long ll;
ll qul_mul(ll a, ll b, ll n) //快速积运算，快速求出a^m^*a^m^并且能防止爆掉
{
ll num = 0;
while (b)
{
if (b&1)
{
num=(num+a)%n;
}
a = (a+a)%n;
b>>=1;
}
return num;
}
ll qul_pow(ll a, ll b, ll n) //快速幂计算快速求出a^m^的值
{
ll sum = 1;
while (b)
{
if (b&1)
{
sum = sum * a % n;
}
a=a*a%n;
b>>=1;
}
return sum;
}
bool Miller_Rabin(ll n)
{
int m = n-1;
int t = 0;
if (n==2)
{
return true;
}
else if (n<2 || !(n&1))
{
return false;
}
//求出n-1 = m*2^t的m和t。
while (!(m&1))
{
m>>=1;
t++;
}
for (int i=0; i<20; i++)
{
//随机数a
ll a = rand()%(n-1) + 1;
// 计算a^m
ll x = qul_pow(a, m, n), y;
for (int j=0; j<t; j++)
{
y = qul_mul(x, x, n); //进行(x*x)%n操作。
//不满足二次探测定理，也就是y得1了但是x并不等于1或者n-1，那么n就一定不是质数。
if (y==1 && x!=1 && x!=n-1)
{
return false;
}
x = y;
}
//上面循环跑完了，如果最后x都不等于1，那么一定是一个合数了。
if (x!=1)
{
return false;
}
}
return true;
}
int main()
{
ll n;
while (cin >> n)
{
cout << Miller_Rabin(n) << endl;
}
return 0;
}
代码
注意在进行素数判断的时候需要用到快速幂。。
#include<cstdio>
#define LL long long
inline int read() {
char c = getchar(); int x = 0, f = 1;
while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return x * f;
}
int N, M, Test[10] = {2, 3, 5, 7, 11, 13, 17};
int pow(int a, int p, int mod) {
int base = 1;
for(; p; p >>= 1, a = (1ll * a * a) % mod)
if(p & 1) base = (1ll * base * a) % mod;
return base % mod;
}
bool Query(int P) {
if(P == 1) return 0;
int t = P - 1, k = 0;
while(!(t & 1)) k++, t >>= 1;
for(int i = 0; i < 4; i++) {
if(P == Test[i]) return 1;
LL a = pow(Test[i], t, P), nxt = a;
for(int j = 1; j <= k; j++) {
nxt = (a * a) % P;
if(nxt == 1 && a != 1 && a != P - 1) return 0;
a = nxt;
}
if(a != 1) return 0;
}
return 1;
}
main() {
N = read(); M = read();
while(M--) puts(Query(read()) ? "Yes" : "No");
}
