Miller-Rabin 素性检验是一种素数判定法则，利用随机化算法判断一个数是合数还是可能是素数。卡内基梅隆大学的计算机系教授Gary Lee Miller首先提出了基于广义黎曼猜想的确定性算法，由于广义黎曼猜想并没有被证明，其后由以色列耶路撒冷希伯来大学的Michael O. Rabin教授作出修改，提出了不依赖于该假设的随机化算法。
算法原理
费马小定理
Miller-Rabin 算法的依据是费马小定理: 假如 $a$ 是一个整数，$p$ 是一个质数，那么 $a^p$ 是 $p$ 的倍数，可以表示为
$$ a^{p-1} \equiv 1(mod \; p) $$
证明过程去 Wiki 上面看吧。
那是不是当一个数 $p$ 满足任意 $a$ ，使得 $a^{p-1} \equiv 1(mod \; p)$ 成立时，它就一定是素数呢？答案是否定的，因为有些合数也满足费马小定理，比如 卢卡斯-卡米切尔数。
虽然说费马小定理有一定概率误判，但如果我们通过一些操作，将这个误判的概率降低至非常小呢？比如对于一个数，我们的算法可以保证 99.99999% 的概率做出正确的判断，那么这个算法不也是很优秀的吗？
所以我们有了 Miller-Rabin 算法。
二次探测定理
Miller-Rabin 算法在费马小定理的基础上加了一个二次探测定理:
若 $p$ 是一个奇素数，则由费马小定理可知: 对于任意 $0 < x < p$ ，有 $x^2 \equiv 1(mod \; p)$ 。由此可推出 $x \equiv \pm 1 (mod \; p)$ 。
证明如下:
$$ x^2 \equiv 1(mod \; p) $$
$$ x^2 - 1 \equiv 0(mod \; p) $$
$$ (x + 1) (x - 1) \equiv 0(mod \; p) $$
由以上式子可知:
$$ (x + 1) \equiv 0(mod \; p) 或 (x - 1) \equiv 0(mod \; p) $$
即
$$ x \equiv \pm 1(mod \; p) $$
二次探测定理的前提是这个数满足费马小定理，所以它可以看作是对费马测试的一个补充，减少它误判概率的。同费马小定理，我们这里用的也是它的逆否命题，即不满足二次探测定理的数一定不是素数。
算法步骤
上面介绍了 Miller-Rabin 算法的两个重要的定理，下面我们来讨论一下它是怎么通过这两个定理来判断一个数是否为素数的。
因为 2 是唯一一个偶数素数，所以它需要单独判断。下面所说的素数都是不包括 2 的。
我们假设需要判断的数是 $p$ 。由于 $p$ 是一个奇数，所以 $p - 1$ 是一个偶数。那么它就可以被分解为 $2^u * r$ 的形式(其中 $u, r$ 都是整数且 $r$ 为奇数)。
为什么要分解成这种形式呢？不急，我们先继续往下看。
根据费马小定理，如果 $p$ 是素数，那么一定有 $a^{2^k * t} \equiv 1 (mod \; p)$ 。这时候是不是很明显了？回头看看刚才的二次探测定理，把 $a^{2^k * t}$ 提出一个平方，不就刚好是它的形式吗？
所以我们随机选择一个数 $a$ ，计算出 $a^t (mod \; p)$ ，然后让其不断自乘，同时结合二次探测定理进行判断。
如果我们自乘后的数模 $p$ 的结果为 1 ，而之前的数模 $p$ 不为 $\pm 1$ ，就表明数 $p$ 是个合数(不满足二次探测定理)
如此循环 $k$ 次，最后得到的数就是 $a^{p-1}$ 。如果这个数模 $p$ 的结果不为 1 的话，就表明 $p$ 是合数(不满足费马小定理)
至此，我们就完成了一次素性探测了。当然，一次探测的误判率还是很高。所以我们需要多判断几次。
正确率
网上找的资料显示: 经过独立 $k$ 轮的 Miller-Rabin 算法误判的概率大约为 $1 / 4^k$ 。
代码实现
素数判断时需要用到快速幂。这个比较简单，不懂的可以看我之前写的快速幂
C++ 实现
由于我是用 C++ 实现的，大整数类懒得写上去，所以使用的时候要注意越界问题。(其实写好大整数类后把数据类型改一下就行了)
Miller-Rabin.h
#pragma once
typedef long long long_t;
inline long_t multiply_mod(const long_t& factor_1, const long_t& factor_2, const long_t& m = 0);
inline long_t quick_pow_mod(const long_t& _base, const long_t& _index, const long_t& m = 0);
bool MillerRabin(const long_t& p, const int& times);Miller-Rabin.cpp
#include "Miller-Rabin.h"
#include <ctime>
#include <cstdlib>
/**
* \brief 模意义下的大整数乘法
* \param factor_1 因数 1
* \param factor_2 因数 2
* \param m 取模, 值为0表示不取模
*/
inline long_t multiply_mod(const long_t& factor_1, const long_t& factor_2, const long_t& m)
{
long_t a = factor_1 % m;
long_t b = factor_2 % m;
long_t res = 0;
// 和快速幂模一样的原理
while (b != 0)
{
if ((b & 1) == 1)
res = (m == 0 ? res + a : (res + a) % m);
a <<= 1;
if (m != 0 && a > m)
a %= m;
b >>= 1;
}
return res;
}
/**
* \brief 一个简单的快速幂模
* \param _base 底数
* \param _index 指数
* \param m 取模, 防溢出, 值为0表示不取模
*/
inline long_t quick_pow_mod(const long_t& _base, const long_t& _index, const long_t& m)
{
long_t index = _index;
long_t base = _base;
long_t res = 1;
while (index != 0)
{
if ((index & 1) == 1)
res = multiply_mod(res, base, m);
index >>= 1;
base = multiply_mod(base, base, m);
}
return (m == 0 ? res : res % m);
}
/**
* \brief Miller-Rabin 素性探测
* \param p 需要判断的数
* \param times 判断次数, 减少误差用
*/
bool MillerRabin(const long_t& p, const int& times)
{
// 判断特殊情况
if (p == 2)
return true;
if (p < 2 || p % 2 == 0)
return false;
long_t r = p - 1;
while (r % 2 == 0)
r /= 2;
// 播种子
srand(static_cast<unsigned>(time(nullptr)));
for (auto i = 0; i < times; i++)
{
long_t a = rand() % (p - 1) + 1;
long_t x = quick_pow_mod(a, r, p);
for (auto u = 1; u * r < p; u *= 2)
{
const long_t tmp = multiply_mod(x, x, p);
// 不满足二次探测定理
if (tmp == 1 && x != 1 && x != p - 1)
return false;
x = tmp;
}
if (x != 1)
return false;
}
return true;
}Python 实现
import random
def MillerRabin(num, times):
if num == 2:
return True
if num < 2 or num % 2 == 0:
return False
r = num - 1
k = 0
while r % 2 == 0:
r /= 2
k += 1
for t in range(times):
a = random.randrange(2, num - 1)
x = pow(a, r, num)
for i in range(1, k + 1):
tmp = x * x % num
if tmp == 1 and x != 1 and x != num - 1:
return False
x = tmp
if x != 1:
return False
return True
def isPrime(num):
if num<2:
return False
lowPrimes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
if num in lowPrimes:
return True
for prime in lowPrimes:
if num % prime==0:
return False
return MillerRabin(num, 8)
if __name__ == '__main__':
while True:
num = input()
print (isPrime(int(num)))Referrer
https://www.cnblogs.com/zwfymqz/p/8150969.html#_labelTop
https://en.wikipedia.org/wiki/Fermat%27s_little_theorem
https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
