[模板] Miller-Rabin 素性测试

背景

在十二省联考中突然忘记 Miller-Rabin 怎么写，决定滚回来写这篇博客

费马小定理

假设 $x$ 是一个整数，$p​$ 是一个质数，那么满足
$$
x^{p-1} \equiv 1 \pmod p
$$
这个定理也是求模质数逆元的依据，即
$$
x^{p-2} \times x \equiv x^{p-1} \equiv 1 \pmod p
$$
因此我们可以随机一些数 $x$ ，判断 $x^{p-1}$ 对 $p$ 取模是否同余于 $1$，如果不是就是一个合数

这种方法已经可以正确判断很多数了，但有一些具有特殊性质的数依旧无法判断

二次探测定理

假设存在一个质数 $p$ 和一个整数 $x​$
$$
x^2 \equiv 1 \pmod p
$$
那么必定满足 $x \equiv 1 \pmod p$ 或 $x \equiv p – 1 \pmod p$

那么我们可以将 $p-1$ 写成 $2^u \cdot v$ 的形式，枚举一个大于 $1$ 小于 $p$ 的正整数 $a$

判断
$$
a^v, a^{2v}, a^{2^2v} \cdots, a^{2^{u-1}v}
$$
中是否有模 $p$ 同余于 $p-1$ 的数，有说明 $p$ 有可能是一个质数，没有就不是一个质数

这个算法的正确率是 $\frac{3}{4}$ ，因此我们可以随机 $k$ 个数用于判定，正确率是 $1 – \left( \frac{1}{4} \right)^k$

当 $k$ 大于等于 $10$ 时，算法错误率小于 $0.00001\%$

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <ctime>
#include <algorithm>

using namespace std;
const int TIMES = 10;

inline long long mul(long long a, long long b, long long mod){
    long long ans = 0;

    for(; b; b >>= 1){
        if(b & 1) ans = (ans + a) % mod;
        a = (a + a) % mod;
    }

    return ans;
}

inline long long powx(long long a, long long b, long long mod){
    long long ans = 1;

    for(; b; b >>= 1){
        if(b & 1) ans = mul(ans, a, mod);
        a = mul(a, a, mod);
    }

    return ans;
}  

inline bool check(long long a, long long x){
    long long tmp = x - 1;
    int k = 0;

    while(!(tmp & 1)){
        tmp >>= 1;
        k++;
    }

    tmp = powx(a, tmp, x);

    if(tmp == 1) return true;

    for(register int i=1; i<=k; i++){
        if(tmp == x - 1) return true;
        tmp = mul(tmp, tmp, x);
    }

    return false;
}


inline bool miller_rabin(long long x){
    if(x < 2) return false;
    if(x == 2 || x == 3 || x == 5 || x == 7) return true;
    if(x % 2 == 0 || x % 3 == 0 || x % 5 == 0 || x % 7 == 0) return false;

    for(register int i=1; i<TIMES; i++){
        long long a = rand() % (x - 2) + 2;
        if(!check(a, x)) return false;
    }

    return true;
}