首先介绍一下miller_rabin算法。
miller_rabin是一种素性测试算法,用来判断一个大数是否是一个质数。
miller_rabin是一种随机算法,它有一定概率出错,设测试次数为
s
s
s,那么出错的概率是
4
−
s
4^{-s}
4−s,至于为什么我也不会证明。我觉得它的复杂度是
O
(
s
l
o
g
2
n
)
O(slog^2n)
O(slog2n),因为你要进行
s
s
s次,每次要进行一次快速幂,每次快速幂要
l
o
g
n
logn
logn次快速乘,每次快速乘又是
l
o
g
n
logn
logn的,所以这一部分是
l
o
g
2
n
log^2n
log2n的,另一部分是把
n
n
n分解成
u
∗
2
k
u*2^k
u∗2k,复杂度是
l
o
g
n
logn
logn,所以k是
l
o
g
n
logn
logn量级的,对于
k
k
k,每次要快速乘,所以这一部分的复杂度也是
l
o
g
2
n
log^2n
log2n的,于是总复杂度
O
(
m
l
o
g
2
n
)
O(mlog^2n)
O(mlog2n)。
下面开始介绍miller_rabin的做法。
首先我们知道,根据费马小定理,如果
p
p
p为质数且
a
a
a与
p
p
p互质,那么有
a
p
−
1
=
1
(
m
o
d
p
)
a^{p-1}=1(mod\ p)
ap−1=1(mod p),如果我们让
a
<
n
a<n
a<n,那么只需要满足
p
p
p为质数即可。但是反过来,满足
a
p
−
1
=
1
(
m
o
d
p
)
a^{p-1}=1(mod\ p)
ap−1=1(mod p)的p不一定是质数。但是幸运的是,对于大多数的
p
p
p,费马小定理的逆定理是对的,但是为了让我们得到正确是结果,我们需要提高正确率。所以接下来我们需要二次探测。
若
n
n
n为质数,并且
x
2
=
1
(
m
o
d
n
)
x^2=1(mod\ n)
x2=1(mod n),那么
x
=
1
x=1
x=1或
x
=
n
−
1
x=n-1
x=n−1,因为
x
2
=
1
x^2=1
x2=1,则
x
=
1
x=1
x=1或
x
=
−
1
x=-1
x=−1,而在模意义下这个
−
1
-1
−1就成了
n
−
1
n-1
n−1。我们来证明一下为什么在
x
2
=
1
(
m
o
d
p
)
x^2=1(mod\ p)
x2=1(mod p),且
x
≠
1
(
m
o
d
p
)
x≠1(mod\ p)
x=1(mod p),且
x
≠
p
−
1
(
m
o
d
p
)
x≠p-1(mod\ p)
x=p−1(mod p)时
p
p
p不是质数。
x
2
=
1
(
m
o
d
p
)
x^2=1(mod\ p)
x2=1(mod p)
x
2
−
1
=
0
(
m
o
d
p
)
x^2-1=0(mod\ p)
x2−1=0(mod p)
(
x
−
1
)
(
x
+
1
)
=
0
(
m
o
d
p
)
(x-1)(x+1)=0(mod\ p)
(x−1)(x+1)=0(mod p)
所以有
p
∣
(
x
−
1
)
p|(x-1)
p∣(x−1)或
p
∣
(
x
+
1
)
p|(x+1)
p∣(x+1)或
p
∣
(
x
+
1
)
(
x
−
1
)
p|(x+1)(x-1)
p∣(x+1)(x−1)
因为
x
≠
p
−
1
(
m
o
d
p
)
x≠p-1(mod\ p)
x=p−1(mod p)
所以
(
x
+
1
)
%
p
≠
0
(x+1)\%p≠0
(x+1)%p=0,
p
p
p不整除
x
+
1
x+1
x+1
因为
x
≠
1
(
m
o
d
p
)
x≠1(mod\ p)
x=1(mod p)
所以
(
x
−
1
)
%
p
≠
0
(x-1)\%p≠0
(x−1)%p=0,
p
p
p不整除
x
−
1
x-1
x−1
所以只能是
p
∣
(
x
+
1
)
(
x
−
1
)
p|(x+1)(x-1)
p∣(x+1)(x−1)
由此我们可以发现
x
+
1
x+1
x+1和
x
−
1
x-1
x−1都不含
p
p
p的所有因子,而它们的乘积却含有的
p
p
p的全套因子,那么在相乘时
x
+
1
x+1
x+1和
x
−
1
x-1
x−1分别提供了一个非
1
1
1非
p
p
p的因子,才能使乘积是
p
p
p的倍数(因为
p
∣
(
x
+
1
)
(
x
−
1
)
p|(x+1)(x-1)
p∣(x+1)(x−1)嘛),当然这两个因子可能相同。
那么,我们可以证明出
p
p
p可以表示为两个非
1
1
1非
p
p
p的数的乘积,那么就证明了
p
p
p是合数而不是质数了。
我们可以根据这个性质进行二次探测。
如果我们设要测试的数为
n
n
n,若
n
n
n是
2
2
2或
n
n
n是偶数,那么我们可以直接判断
n
n
n的奇偶,否则令
p
=
n
−
1
p=n-1
p=n−1,将
p
p
p分解成
u
∗
2
k
u*2^k
u∗2k,枚举
k
k
k来不断进行二次探测。则那么我们随机一个底数
x
x
x,快速幂求出
x
u
(
m
o
d
n
)
x^u(mod\ n)
xu(mod n),然后不断的
x
∗
x
m
o
d
p
x*x\ mod\ p
x∗x mod p来进行二次探测。
至于为什么在代码中这么做是对的,我看到以了个很好的解释,在这里分享一下。
a
p
−
1
=
1
(
m
o
d
p
)
a^{p-1}=1(mod\ p)
ap−1=1(mod p)
a
u
∗
2
k
=
1
(
m
o
d
p
)
a^{u*2^k}=1(mod\ p)
au∗2k=1(mod p)
(
a
u
∗
2
k
−
1
)
∗
(
a
u
∗
2
k
−
1
)
=
1
(
m
o
d
p
)
(a^{u*2^{k-1}})*(a^{u*2^{k-1}})=1(mod\ p)
(au∗2k−1)∗(au∗2k−1)=1(mod p)
那么我们把
a
u
∗
2
k
−
1
a^{u*2^{k-1}}
au∗2k−1看作
x
x
x,那么就是在当前
a
u
∗
2
k
=
1
(
m
o
d
p
)
a^{u*2^k}=1(mod\ p)
au∗2k=1(mod p)的时候验证
a
u
∗
2
k
−
1
(
m
o
d
p
)
a^{u*2^{k-1}}(mod\ p)
au∗2k−1(mod p)是否是
p
−
1
p-1
p−1或者
1
1
1即可。
另外最后别忘了用费马小定理来判断一下,代码中的
x
x
x其实的
x
p
−
1
(
m
o
d
p
)
x^{p-1}(mod\ p)
xp−1(mod p)之后的结果,所以只需要看最后的
x
x
x是否是
1
1
1即可。
最后是代码:
#include <bits/stdc++.h>
using namespace std;
int q,m,s=5;//s为测试次数(选了几次底数a)
inline long long ksc(long long x,long long y,long long mod)//快速乘
{
long long res=0;//注意赋初值
// x%=mod;
while(y)
{
if(y&1)
res=(res+x)%mod;
x=(x<<1)%mod;
y>>=1;
}
return res;
}
inline long long ksm(long long x,long long y,long long mod)
{
long long res=1;//注意设为1,不是0
// x%=mod;
while(y)
{
if(y&1)
res=ksc(res,x,mod);
x=ksc(x,x,mod);
y>>=1;
}
return res;
}
inline int miller_rabin(long long n)
{
if(n==2||n==3||n==5||n==7||n==11)
return 1;
if(n<2||!(n%2)||!(n%3)||!(n%5)||!(n%7)||!(n%11))
return 0;
long long x,pre,u;//pre为上次的结果
int k=0;//k为n分解成了2的多少次方
//最终n被分解为u*2^k
u=n-1;//求x^u%n;
while(!(u&1))//u为偶数则右移,否则就停
{
++k;
u>>=1;
}
srand(time(0)+19260817);
for(int i=1;i<=s;++i)
{
x=rand()%(n-2)+2;//生成一个[2,n)的随机底数
x=ksm(x,u,n);//先求出x^u mod n
pre=x;
for(int j=1;j<=k;++j)//把移位减掉的量补上,并在这地方进行二次探测
{
x=ksc(x,x,n);
if(x==1&&pre!=1&&pre!=n-1)//二次探测定理,这里如果x = 1则pre 必须等于 1,或者 n-1,否则可以判断n不是素数
return 0;
pre=x;
}
if(x!=1)//费马小定理
return 0;
}
return 1;
}
int main()
{
scanf("%d%d",&q,&m);
for(int i=1;i<=m;++i)
{
int x,pd;
scanf("%d",&x);
pd=miller_rabin(x);
if(pd==1)
printf("Yes\n");
else
printf("No\n");
}
return 0;
}
最后发一个题单:
洛谷3383
poj1811
HDU2138
(都是模板题)
最后:特别鸣谢wpc、was_n、DT_Kang和liuzhangfeiabc帮忙给出的一些证明