Miiler Rabin 算法

Miller Rabin 算法

参考：【朝夕的ACM笔记】数论-Miller Rabin素数判定 - ~(o°ω°o)

作用在于判断素数，它有前置两个定理

费马小定理

Theorem 设 \(p\) 是一个素数，\(a\in \Z^+\) 且不是 \(p\) 的倍数，那么有 \(a^{p-1}\equiv1(\bmod p)\)

Proof 我们知道，若有 \(i\) 和 \(j\) 在模 \(p\) 意义下不同余，那么 \(i\times a\) 和 \(j\times a\) 在模 \(p\) 意义下也不同余

那么有： \(1\times 2\times 3\times\dots\times(p-1)\equiv(1\times a)\times(2\times a)\times\dots\times(p-1)*a(\bmod p)\)

设 \(W=\Pi_{i=1}^{p-1}\) ，那么 \(W\equiv W\times a^{p-1}(\bmod p)\)

又知：\(\gcd(W,p)=1\), 所以，\(a^{p-1}\equiv 1\pmod p\)

证毕。

我们可以猜想，费马小定理的逆定理是否成立

也就是说，如果 \(a\in \Z^+\) 且不是 \(p\) 的倍数，而且 \(a^{p-1}\equiv1\pmod p\) ，是否能说 \(p\) 是一个素数？

答案是不能，反例：对 \(a = 2\)，满足 \(2^{n−1} \bmod n = 1\) 的非素数 \(n\) 是存在的，比如 \(n = 341 = 11 × 31\)

伪素数：对于整数 \(a\)，称满足 \(a^{n−1} \bmod n = 1\) 的合数为以 \(a\) 为底的伪素数。

经测试，前 10 亿的自然数中，同时以 2 和 3 为底的伪素数有 1272 个。

我们用费马小定理验证素数的话，出错的概率大概只有 0.000025。

于是我们初步得到了一个判断素数的算法：

Solution 随机选取若干个小于待测整数的正整数作为底 a，然后用费马小定理来测试

但是，有一些坑爹的合数能通过所有的测试，也就是 \(\text{Carmichael}\) 数

二次探测定理

Theorem 若 \(p\) 为素数, $a^2≡1\pmod p $，那么 \(a≡±1\pmod p\)

Proof 其实感觉是显然的。由题，\(p|(a^2-1)\) ，即 \(p|(a-1)(a+1)\)

那么，因为 \(p\) 是素数，根据唯一分解定理，\(p|(a-1)\) 或者 \(p|(a+1)\) 。得证。

所以，如果 $a^2≡1\pmod p $ ，但是不满足 \(a≡±1\pmod p\) ，那么 \(p\) 不是素数。

MR

二次探测定理和素数判定有什么用呢？

首先，根据 \(Miller Rabin\) 算法的过程

假设需要判断的数是 \(p\) （假设 \(p\) 是奇数，是偶数时显然不是素数）

我们把 \(p−1\) 分解为 \(2^k+t\) 的形式（当 \(p\) 是素数，\(a^{2^k+t}\equiv 1\pmod p\)）

我们随机一个 \(a\) ，计算出 \(a^t\) ，让它不断的自乘，

如果在某时，\(a^{2k}\equiv1\mod p\) ，但是不满足 \(a^k≡±1\pmod p\) ，那么 \(p\) 不是素数

如果最后算到 \(a^{2^k+t}\bmod p\ne1\) 那 \(p\) 不是素数。

• 正确性：首先是底数的选取：Miller Rabin检测依旧有错误的概率，虽然微乎其微，但我们显然不想接受。但通过选取适合的底数，可以避免这一情况的发生。

  在 \(int\) 范围内，选取 \(2,7,61\) 三个数作为底数可以保证 \(100\%\) 正确。

  在 \(long long\) 范围内，选取 \(2,325,9375,28178,450775,9780504,1795265022\) 七个数作为底数可保证 \(100\%\) 正确。

  故正常情况下时间复杂度最多为 \(O(n \log^3 n)\) ，常数为 \(7\)。

代码实现：

#include<cstdio>
#define LL long long 
inline int read() {
    char c = getchar(); int x = 0, f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * f;
}
int N, M, Test[10] = {2, 3, 5, 7, 11, 13, 17};
int pow(int a, int p, int mod) {
    int base = 1;
    for(; p; p >>= 1, a = (1ll * a * a) % mod) 
        if(p & 1) base = (1ll * base * a) % mod;
    return base % mod;
}
bool Query(int P) {
    if(P == 1) return 0;
    int t = P - 1, k = 0;
    while(!(t & 1)) k++, t >>= 1;
    for(int i = 0; i < 4; i++) {
        if(P == Test[i]) return 1;
        LL a = pow(Test[i], t, P), nxt = a;
        for(int j = 1; j <= k; j++) {
            nxt = (a * a) % P;
            if(nxt == 1 && a != 1 && a != P - 1) return 0;
            a = nxt;
        }
        if(a != 1) return 0;
    }
    return 1;
}
main() { 
    N = read(); M = read();    
    while(M--) puts(Query(read()) ? "Yes" : "No");
}