Miller Rabin算法学习笔记

定义:##

Miller Rabin算法是一个随机化素数测试算法,作用是判断一个数是否是素数,且只要你脸不黑以及常数不要巨大一般来讲都比(O(sqrt n))的朴素做法更快。

定理:##

Miller Rabin主要基于费马小定理:

[a ^ {p-1} equiv 1 (mod p)$$其中$p$是质数。 于是就有~~闲得没事干的~~一群科学家们想,这个问题的逆命题是否成立呢? > 逆命题:若对于任意$a$,$a ^ {p-1} equiv 1 (mod p)$都成立,那么$p$是质数。 在很长一段时间里,所有人几乎都以为它是成立的。~~然鹅你们手玩一个$a=8, p = 9$试试~~ 是的,这个东西被搞出了反例。不过幸运的是,用这个办法测试通过的数,还是有很大概率是质数的。 这好办,我们多搞几次不就可以当做它就是质数了吗!~~脸黑另说~~ ##算法流程:## 首先我们还得了解一个叫二次探测定理的东西: > $$若p是质数,且x^2 equiv 1 (mod p), 则有x equiv ±1 (mod p)]

证明很简单,第一个式子右边丢过去平方差即可。由于p是质数,所以它肯定不是((x-1)和(x+1))凑起来的,故两个里面总有一个是(p)的倍数。
而且很容易脑补的是,这个东西的逆命题是成立的。(划重点)
所以根据这两个定理,我们设计一波算法:
假设我们要判断的数是(p),那么(2)特判一波,剩下的质数肯定是奇数。
所以(p-1)一定是一个偶数。然后就好办啦!
我们把(p-1)分解成(2^k * t),当(p)是素数时,根据费马小定理有$$a ^ {2^k * t} equiv 1 (mod p)$$
那么我们随机出一个(a),然后求出(a^t),再不断乘上(a),每次进行二次探测,边乘边模,若乘之前不符合二次探测,而乘之后符合,那么p是合数,不符合题意。自乘(k)次,最后得到(a^{p-1}),如果模(p)不等于1,则也是合数。(不符合费马小定理)
老祖宗告诉我们(这个我也不会证),每一次通过测试的数不是质数的概率为(frac{1}{4}),则测试(k)次,错误的概率为(frac{1}{4^k})(k>6)的时候基本就血赚了。

代码:##

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int c[23] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43};
int n, m;
inline ll read() {
    ll cnt = 0, f = 1; char c;
    c = getchar();
    while (!isdigit(c)) {if (c == '-') f = -f; c = getchar();}
    while (isdigit(c)) {cnt = (cnt << 3) + (cnt << 1) + c - '0'; c = getchar();}
    return cnt * f;
}
inline ll ksm(ll a, ll b, ll c) {
    ll ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % c; 
        a = a * a % c, b >>= 1;
    }
    return ans % c;
}
bool miller_rabin(int p) {
    if (p == 1) return false;
    if (p == 2) return true;
    if (p % 2 == 0) return false;
    bool f = 1;
    for (register int i = 0; i <= 13; ++i) {
        if (c[i] == p) return true;
        ll x = p - 1, y = 0;
        while (x % 2 == 0) x /= 2, ++ y;   // 将p-1分解成2^y*x
        ll cur = ksm(c[i], x, p);  //计算出a^x % p
        if (cur == 1) continue; //小优化,如果此时结果为1,那么无论如何自乘也为1 
        for (register int j = 1; j <= y; ++j) {
            ll nxt = cur * cur % p;  //不断自乘 
            if (nxt == 1 && cur != p - 1 && cur != 1) {
                f = 0;
                break;
            }
            cur = nxt;
        }
        if (cur != 1) f = 0;
        if (!f) break;
    }
    return f;
}
int main() {
	n = read(); m = read();
	while (m--) {printf(miller_rabin(read()) ? "Yes
" : "No
");}
    return 0;
}


原文地址:https://www.cnblogs.com/kma093/p/11115535.html