Miller-rabin 素性测试
无名之雾
·
·
算法·理论
前置
费马小定理
设 n 是素数，a 是正整数且与 n 互素，那么有 a^{n-1} \equiv 1(\bmod n)。
费马素性测试
为了测试 n 是否为素数，我们可以在 1 与 n 之间任选一个随机的基数 a。
-
如果 a^{n-1} \equiv 1(\bmod n) 不成立，那么根据费马小定理的逆命题，n 肯定不是素数。
-
如果 a^{n-1} \equiv 1(\bmod n) 成立，那么 n 大概率 是素数。尝试的 a 越多，正确的概率越大。特别的，有一些合数对任何 a 都成立，这种数被称为 Carmichael 数。
二次探测定理
如果 p 是一个奇素数，且 e \ge 1，则方程 x^2 \equiv 1(\bmod p^e) 仅有两个解，x=1 和 x=-1。当 e=1 时，方程 x^2 \equiv 1(\bmod p) 仅有两个解，x=1 和 x=p-1。
我们称 x=1 和 x=p-1 为 x 对模 p 来说 1 的平凡平方根。
证明
## Miller-rabin 素性测试
Miller-rabin 用到这个方程，$x^2 \equiv 1(\bmod n)$。如果一个数 $x$ 满足上述方程，但 $x$ 不等于平凡平方根 $1$ 或 $n-1$。那么称 $x$ 是对模 $n$ 来说 $1$ 的非平凡平方根。
根据二次探测定理的逆命题，如果存在 $x$ 是对模 $n$ 来说 $1$ 的非平凡平方根，则 $n$ 是合数。
### Miller-Rabin 素性测试的步骤
输人 $n>2$，且 $n$ 是奇数,测试它是否为素数。
根据费马测试,如果 $a^{n-1} \equiv 1(\bmod n)$ 不成立,那么 $n$ 肯定不是素数。
由于 $n-1$ 太大,显然不能直接计算 $a^{n-1}$。所以需要用到一个数学技巧。把 $n-1$ 表示为幂的形式,然后借助快速冠来计算。
令 $n-1=2^tu$，其中 $u$ 是奇数,是正整数。求 $t$ 和$u$ 编程时这样做。$n-1$ 的二进制表示是奇数 $u$ 的二进制表示后面加个 $0$。
选一个随机的基值 $a$，有
$$a^{n-1} \equiv (a^u)^{2^t} (\bmod n)$$
下面就可以用快速幂来计算了。为了计算 $a^{n-1}(\bmod n)$，可以先算出 $a^{u}(\bmod n)$，然后对结果连续平方 $t$ 次取模。这是因为符合乘法模运算规则 $(c\times d)\bmod n=(c\bmod n\times d \bmod n)\bmod n$。
### Miller-Rabin 素性测试的时间复杂度和出错概率
Miller-Rabin 素性测试需要用多个随机的基值 $a$ 来做以上的测试。设有 $s$ 个 $a$，共做 $s$ 次测试，出错的概率为$2^{-s}$。当 $s=35$ 时,出错概率已经小到可以忽略不计了。
计算复杂度分析，算法做了 $s$ 次运算,再考虑到 $t$ 次快速幂，总复杂度在最坏情况下为 $O(s\times(\log n^3))$。
## 代码
```cpp
#include<bits/stdc++.h>
using namespace std;
#define int long long
inline int read(){
int s=0,w=1;
char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
while(ch>='0'&&ch<='9') s=s*10+ch-'0',ch=getchar();
return s*w;
}
inline void write(int x){
if(x==0){putchar('0');return;}
int len=0,k1=x,c[10005];
if(k1<0)k1=-k1,putchar('-');
while(k1)c[len++]=k1%10+'0',k1/=10;
while(len--)putchar(c[len]);
}
int fpow(int x,int y,int m){
int res=1;
x%=m;
while(y){
if(y&1)res=(res*x)%m;
x=(x*x)%m;
y>>=1;
}
return res;
}
bool witness(int a,int n){
int u=n-1,t=0;
while(u&1==0)u>>=1,t++;
int x1=fpow(a,u,n);
for(int i=1;i<=t;i++){
int x2=fpow(x1,2,n);
if(x2==1&&x1!=1&&x1!=n-1)return 1;
x1=x2;
}
if(x1!=1)return 1;
return 0;
}
bool miller_rabin(int n,int s){
if(n<2)return 0;
if(n==2)return 1;
if(n%2==0)return 0;
for(int i=0;i<s&&i<n;i++){
int a=rand()%(n-1)+1;
if(witness(a,n))return 0;
}
return 1;
}
signed main(){
srand(time(0));
int l=read(),r=read(),cnt=0;
for(int i=l;i<=r;i++){
if(miller_rabin(i,35))cnt++;
}
write(cnt);
return 0;
}
```