一、算法简介
在算法竞赛中,我们时常会遇到需要判断一个数是否为质数的问题。我们常常利用筛法来解决这个问题,但是当需要判断的数变得很大时,筛法已经无法满足我们的需求。于是我们采用了一个新的方法:Miller-Rabin素数检测。
二、算法分析
1.前置知识
(1)费马小定理
- 由费马小定理可知,若p pp为质数且a aa不是p pp的倍数,a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(mod\ p)ap−1≡1(mod p)成立。
- 反过来,如果存在某个p pp满足:a p − 1 ≡ 1 ( m o d p ) a^{p-1}\equiv1(mod\ p)ap−1≡1(mod p),并不一定意味着p pp是质数,例如2 341 − 1 ≡ 1 ( m o d p ) 2^{341-1}\equiv1(mod\ p)2341−1≡1(mod p) ,但341 341341是合数。
- 满足该同余等式的合数被称为费马伪素数。但这样的数并不多,满足这个条件的数很多都是质数。
(2)另一个定理
①定理内容
- 如果p pp是一个素数,0 < x < p 0<x<p0<x<p, 则方程 x 2 ≡ 1 ( m o d p ) x^{2}\equiv1(mod\ p)x2≡1(mod p)的解为x = 1 x=1x=1或x = p − 1 x=p-1x=p−1
- 反之,如果x 2 ≡ 1 ( m o d p ) x^2\equiv1(mod\ p)x2≡1(mod p)的解不是x = 1 x=1x=1或x = p − 1 x=p-1x=p−1 那p pp就不是素数
②定理证明
x 2 ≡ 1 ( m o d p ) x^2\equiv1(mod\ p)x2≡1(mod p) 可以移项变成x 2 − 1 ≡ 0 ( m o d p ) x^2-1\equiv0(mod\ p)x2−1≡0(mod p)
展开可得( x − 1 ) ( x + 1 ) ≡ 0 ( m o d p ) (x-1)(x+1)\equiv0(mod\ p)(x−1)(x+1)≡0(mod p),即p ∣ ( x − 1 ) ( x + 1 ) p|(x-1)(x+1)p∣(x−1)(x+1)或x = p − 1 x=p-1x=p−1
由欧几里得引理可得:p pp素数,如果p ∣ b c p|bcp∣bc那么p ∣ b p|bp∣b或者p ∣ c p|cp∣c。所以若p pp是个质数,则p ∣ x − 1 p|x-1p∣x−1或p ∣ x + 1 p|x+1p∣x+1,即x ≡ ± 1 ( m o d p ) x\equiv±1(mod\ p)x≡±1(mod p)
所以解为x = 1 x=1x=1或x = p − 1 x=p-1x=p−1
2.算法过程
假设n nn是一个大于2 22素数,则n − 1 n-1n−1为偶数,可被表示为2 s d 2^sd2sd,s ss和d dd都是正整数且d dd为奇数。 这样一串数字(每一个都是前一个的平方)的性质。
考虑a d , a 2 d , a 4 d , . . . , a 2 s d a^{d},a^{2d},a^{4d},...,a^{2^sd}ad,a2d,a4d,...,a2sd这一串数。我们已经知道对于奇素数n nn,a n − 1 ≡ 1 ( m o d n ) a^{n-1}\equiv1(mod\ n)an−1≡1(mod n)即a 2 s d ≡ 1 ( m o d n ) a^{2^sd}\equiv1(mod\ n)a2sd≡1(mod n)(a aa是x xx倍数的情形特判)。根据第二个定理的内容:
- 若第一个数是取余结果是1 11,则显然后续的所有数都是1 11且取余结果是1 11,满足条件。
- 若除开头和尾,中间出现取余结果为1 11的数,则说明解不是x = 1 x=1x=1或者x = p − 1 x=p-1x=p−1,即这个数不是素数。
- 若非最后一个值中出现取余结果为p − 1 p-1p−1的情况,则说明接下来的所有数都是1 11且取余结果都是1 11,满足条件。
也就是说这串数要么全是1 11,要么中间某个数是− 1 -1−1,并且从那往后全是1 11 。如果满足这个条件,x xx就很可能是素数,否则一定不是质数。于是我们选择一些a aa并验证是否符合条件即可,这就是米勒-拉宾素性检验。
时间复杂度为O ( k log n ) O(k\log{n})O(klogn)
三、模板程序
米勒-拉宾素性检验是一种概率算法,有一定的概率出错。但是,Jim Sinclair发现了一组数:2 , 325 , 9375 , 28178 , 450775 , 9780504 , 1795265022 2,325,9375, 28178,450775,9780504,17952650222,325,9375,28178,450775,9780504,1795265022。用它们做a aa,2 64 2^{64}264以内不会出错,我们使用这组数,就不用担心运气太差了。因为使用固定的7个数,所以时间复杂度为O ( 7 log n ) O(7\log{n})O(7logn),远远优于传统方法的O ( n ) O(\sqrt{n})O(n)。
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <cstdio>
using namespace std;
typedef long long ll;
ll qpow(ll a,ll n,ll p)
{
ll ans=1;
while(n)
{
if(n&1)
ans=(__int128)ans*a%p;
a=(__int128)a * a % p;
n>>=1;
}
return ans;
}
bool is_prime(ll x)
{
if(x<3)
return x==2;
if(x%2==0)
return false;
ll A[7]={2,325,9375,28178,450775,9780504,1795265022},d=x-1,r=0;
while(d%2==0)
d/=2,++r;
for(ll i=0;i<=6;i++)
{
ll v=qpow(A[i],d,x);
if(v<=1 || v==x-1)
continue;
for(ll i=0;i<r;i++)
{
v=(__int128)v*v%x;
if(v==x-1 && i!=r-1)
{
v=1;
break;
}
if(v==1)
return false;
}
if(v!=1)
return false;
}
return true;
}
int main()
{
while(1)
{
ll temp;
scanf("%lld",&temp);
cout<<is_prime(temp)<<endl;
}
return 0;
}