算法思想
对于大于2的素数n，将n-1拆分为
对于待测数n，当找到的a不符合上述条件，那么称a为一个"witness for the compositeness of n"，且n一定是一个合数。否则称a为一个"strong liar"，且n是一个基于a的"strong probable prime"。
准确性
所有的奇合数都有很多的a满足"witness"的条件，不过目前为止还没有确定的算法能够直接根据n生成这样的数a，于是我们可以多次随机抽取1~n-1中的整数并做测试。
当我们k次随机选取a测试时，一个合数被该算法判定为素数的概率是4^(-k)。
一种实现
#include <iostream>
#include <cstdlib>
using namespace std;
typedef long long ll;
ll mod_pow(ll x, ll y, ll m) {
ll base = x, res = 1;
while (y) {
if (y&1) (res*=base)%=m;
(base*=base)%=m;
y>>=1;
}
return res;
}
bool MillerRabin(ll n, int k) {
if (n==2||n==3||n==5||n==7||n==11||n==13) return true;
if (n==1||n%2==0||n%3==0||n%5==0||n%7==0||n%11==0||n%13==0) return false;
ll d=n-1;
int r=0;
while (d%2==0) {
d>>=1;
++r;
}
for (int i=0;i<k;++i) {
ll a=rand()%(n-2)+2;
ll x=mod_pow(a,d,n);
bool flag = false;
if (x==1||x==n-1) continue;
for (int j=0;j<r-1;++j) {
x=mod_pow(x,2,n);
if (x==1) return false;
if (x==n-1) {
flag=true;
break;
}
}
if (flag) continue;
return false;
}
return true;
}
int main() {
ll n;
while (cin>>n) {
if (MillerRabin(n,5))
cout<<"n is a prime number\n";
else cout<<"n is not a prime number\n";
}
return 0;
}
其中ll mod_pow(ll x, ll y, ll m)实现的是模m意义下的快速幂运算。
扩展阅读
- 详细的准确性、时间复杂度分析和证明请参照https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
