Miller-Rabin随机性素数测试方法 [CodeVS 1702] 素数判定2_综合

Miller-Rabin算法用于判定某数x是否为素数。如果x被判定为合数，它一定是合数。如果x被判定为素数，它有很大的概率是素数，此概率取决于参数。

费马小定理：
如果n为素数，那么对于任意与n互质的a，有

a n ? 1 \equiv 1 (mod n) (?)

$a^{n-1} \equiv 1 \pmod n \ \ \ \ \ (*)$

如果a与n不互质，上式一定不成立。因为它意味着 $aa^{n-2} \equiv 1 \pmod n$ ，而a存在模n的乘法逆元当且仅当a, n互质。

如果对于某个a，合数n满足(*)式，那么称n为一个基为a的伪素数。如341是基为2的伪素数。

如果某个合数n对于任意与n互质的a都有(*)式成立，那么称n为Carmichael数。前3个Carmichael数是561, 1105和1729。我们也许会想随机几个a，如果碰到一个与n不互质的基数，就能正确检验了，然而如果n只有大素数因子，这样的基数a很难找到。

只有22个小于10000的基为2的伪素数，只有255个小于100000000的Carmichael数。直接上费马小定理，在实际应用中是可行的，在算法竞赛等领域是不可行的。

Miller-Rabin算法有效地降低了出错的概率。这个概率取决于参数，与被测试的数无关。

设n是待测试的大于1的整数，决定精确程度的参数为s，Miller-Rabin算法的伪代码如下：

Miller-Rabin(n, s)特判n为偶数的情形计算t, u使得n-1=2^t*u且u是奇数，即，t是n-1的二进制表示中末尾0的个数，u是去掉这些0后得到的数for i=1 to s令a为(1, n)内的随机整数x = a^u mod nfor j=1 to t and x != 1y = x*x mod nif y == 1 and x != n-1 // x是以n为模1的非平凡平方根 => n是合数return Compositey = xif x != 1 // a^(n-1) != 1 => n是合数return Compositereturn Prime // 通过素性测试 => n有至少(1-2^(-s))的概率为素数

框架来自《算法导论》，我在不改变算法行为的前提下做了一点改进。

关于算法原理的两点说明：
1) p是奇素数=> $x^2 \equiv 1 \pmod p$ 仅有两平凡解 $x \equiv \pm 1$
证明：方程等价于 $(x+1)(x-1) \equiv 0 \pmod p$ ，由于p是素数，必有p整除(x+1)或(x-1)。p整除(x+1) <=> $x \equiv -1$ ，p整除(x-1) <=> $x \equiv 1$ 。由于p是大于2的奇数， $1 \not \equiv -1$ 。证毕。
由此，当n>1， $x^2 \equiv 1 \pmod n$ 有非平凡解=>n不是奇素数。又，当n=2时，仅有一个平凡解 $x \equiv 1$ ，所以n也不是偶素数。n是合数。

2) 可以证明：在运用原理1)的条件下，如果n是一个奇合数，则n为合数的证据数目至少为(n-1)/2。证据是指判定n为合数的基底a。于是，在[1, n)共(n-1)个数内随机挑选一个，得到证据的概率至少为1/2。把合数误判为素数意味着s次循环都能没找到证据，概率为2^(-s)。

一本由林厚从主编、王宏主审、东南大学出版社出版的《信息学奥赛之数学一本通》告诉我们，最好只生成奇数，以省去原理1)的使用，“进一步提高测试速度”，这当然是错误的。

关于算法实现的两点说明：
1) 计算x = a^u mod n可使用快速幂，整个算法的时间复杂度为 $O(s\lg n)$ 。两个long long类型的数相乘并取模，尽管得到一个long long类型的数，但在相乘的那步就溢出了。解决方案一是采用快速加/快速乘/快速模/慢速模，它们是同一个东西，你懂的。解决方案二是像骆可强《论程序底层优化的一些方法与技巧》一文中提到的那样，借助于浮点运算和最终答案在long long范围内的事实（其实我也没完全理解，暂且当作一个trick记住吧），特判负数是由于对浮点运算四舍五入。

2) 《算法导论》中说，对几乎所有可以想象到的应用，选取s=50应该是足够的。

CodeVS 1702 素数判定2

#include <cstdio>
#include <cstdlib>
#include <ctime>
using namespace std;
typedef long long ll;
const int s = 10;inline ll mul_mod(ll a, ll b, ll n)
{ll c = a*b-(ll)((long double)a*b/n+0.5)*n;return c<0 ? c+n : c;
}ll fast_exp(ll a, ll x, ll n)
{ll b = 1;while (x) {if (x & 1)b = mul_mod(b, a, n);a = mul_mod(a, a, n);x >>= 1;}return b;
}bool MR(ll n)
{if (!(n&1))return n == 2;ll t = 0, u;for (u = n-1; !(u&1); u >>= 1)++t;for (int i = 0; i < s; ++i) {ll a = rand()%(n-2)+2, x = fast_exp(a, u, n);for (ll j = 0, y; x != 1 && j < t; ++j, x = y) {y = mul_mod(x, x, n);if (y == 1 && x != n-1)return false;}if (x != 1)return false;}return true;
}int main()
{srand(time(0));ll p;scanf("%lld", &p);puts(MR(p) ? "Yes" : "No");return 0;
}