米勒拉宾算法——素性测试
背景
对于一个数n，如果想要判断它是否为素数，常规的方法为试除法。即，让n依次除以2到
n
\sqrt{n}
n以内的整数。如果有出现除尽的情况，则为合数。
该方法的时间复杂度为
O
(
n
)
O(\sqrt{n})
O(n)，在面对n为长整型的时候有可能超出时间要求。
因为时间复杂度的优势，对于大整数的素性判定，要使用米勒拉宾算法。
伪素数判定
在此之前介绍一种伪素数判定方法——小费马定理。伪素数判定的意思是，对于大部分的数，可以判断正确，但是对于一些特殊的数字将会判断错误。
小费马定理为：若有素数p，则对任意的数a( a为正整数 ，且a < p)，
a
p
−
1
≡
1
(
m
o
d
p
)
a^{p-1}≡1 (mod \ p)
ap−1≡1(mod p)。反之，若有任意的a( a为正整数 ，且a < p)使得p不满足
a
p
−
1
≡
1
(
m
o
d
p
)
a^{p-1}≡1 (mod \quad p)
ap−1≡1(modp)，p一定为合数。 可以发现若是能够举出所有的a，都能满足上式，是不是就说明p是素数呢？其实不是因为有一类合数也可以做到这一点，这一类合数叫做Carmichael数。前三个这样的数是561, 1105, 1729。这样的数真让人不爽。所以采用这种方法测出来的所谓素数是不一定的。叫做伪素数。哪怕你枚举出所有的a，也不可避免。
更可靠的伪素数测试方法
在小费马定理的基础上有人设计出米勒拉宾随机素数测试法。可以大大的提高检测素数的正确性，但是同样并非一定正确，错误可能性却小到可以接受。
该方法同样利用了枚举多个a的做法，以提高算法的可靠性，对于每一个a，又采用了特殊的方法处理。这基于另外一个定理：如果p是素数，x是小于p的正整数，且
x
2
m
o
d
p
=
1
x^2 \ mod \ p=1
x2 mod p=1，那么要么x=1，要么x=p-1。该定理证明如下：如果p为素数，x是小于p的正整数， 且
x
2
m
o
d
p
=
1
x^2 \ mod \ p=1
x2 mod p=1 ，说明p能够整除**(x+1)和(x-1)。但是p是素数，那么只可能是x-1能被p整除(此时x=1**)或x+1能被p整除(此时 x=p-1)。
多次判定
判断一个数是不是素数光靠上面的方法是不可靠的，因为p如果是合数的话，也有可能有
x
2
≡
1
m
o
d
(
p
)
x^{2}≡1 \ mod (p)
x2≡1 mod(p)且 x=1或者 x =p-1；但是多排除几次p不为合数的话，就增大了p是素数的可能性 ，这是这个算法的核心思想。因此例如341这个数。可知
(
2
340
)
≡
1
(
m
o
d
341
)
( 2^{340} ) ≡ 1 (mod\ 341 )
(2340)≡1(mod 341);
(
2
170
)
≡
1
(
m
o
d
341
)
( 2^{170} ) ≡ 1 (mod\ 341 )
(2170)≡1(mod 341)； 但是发现
2
85
m
o
d
341
≡
32
2^{85} \ mod \ 341 ≡ 32
285 mod 341≡32。这足以证明341是一个合数，而不是一个素数。
首先要判断数n是不是2，再判断n是不是奇数。然后尽可能的在令d=n-1，在d中除去2，使得
n
=
d
(
2
t
)
n=d(2^t)
n=d(2t),d为奇数，t的值并不关心。如果n是一个素数，那么或者
a
d
m
o
d
n
=
1
a^d \ mod \ n=1
ad mod n=1，或者存在某个i使得
a
d
2
i
m
o
d
n
=
n
−
1
,
(
0
≤
i
<
r
)
a^{d2^i} \ mod \ n=n-1,(0\leq i<r)
ad2i mod n=n−1,(0≤i<r)（注意i可以等于0，这就把a^d mod n=n-1的情况统一到后面去了）。
求
a
b
m
o
d
(
n
)
a^b \ mod(n)
ab mod(n)的算法以及求
d
2
d^2
d2的算法是采用的快速幂取模算法。但是在d为long long的情况下有可能乘法溢出。有其他算法可以处理，这里不是重点。
C代码实现
#include <iostream>
#include <cstdio>
using namespace std;
typedef long long LL;
LL mulmod( LL a, LL b , LL p )
{
LL d =1;
a = a%p;
while( b>0 )
{
if(b&1)
d = (d*a)%p;
a = (a*a)%p;
b>>=1;
}
return d;
}
bool witness( LL a,LL n)
{
LL d = n-1 ;
if( n ==2 ) return true ;
if( !(n&1) ) return false ;
while(!(d&1)) d = d/2;
LL t = mulmod(a,d,n);
while((d!=n-1) && (t!=1)&&(t!=n-1))
{
t = mulmod( t ,2,n);
d=d<<1;
}
return (t==n-1)||(d&1);
}
bool isprime( LL n)
{
int a[3] = {2,7,61};
for(int i=0;i<3;i++)
if(!witness(a[i],n))
return false;
return true;
}
int main()
{
LL s;
cin>>s;
if(isprime(s))
cout<<"YES";
else
cout<<"NO";
return 0;
}
PS：入门蒟蒻，如有错误请指正。
多年后的PS: 大学时的一篇水文，没想到还有人看，再次编辑下公式和排版。(2023/11/20)