Miller-Rabin素性测试与Pollard’s rho算法
Miller-Rabin素性测试（Miller–Rabin primality test）
我们知道，朴素的素性测试的实质是枚举因子，我们枚举2至\sqrt{n}的所有正整数，一一检验是否为因子，这样做的复杂度是O(\sqrt{n})。能不能使复杂度更低？
费马小定理
关于费马小定理的阐述参见：数学笔记：数论基础（部分） | KSkun’s Blog。
我们知道了费马小定理的式子，那么它的逆定理如果成立，我们就得到了一种O(\log n)判定方法。可惜的是，逆定理并不成立。
我们选取若干a，每次检查a^{p-1} \bmod p是否为1，这种测试方法叫做费马测试。然而，如果有合数使得选取某个a时通过测试，则称该合数为以a为底的伪素数。能通过所有小于p的底数的合数p称为Carmichael数。如果这些数很大，大过了我们需要的范围，那还好。可是561就是一个很小的Carmichael数，这种测试仍然不能满足需求，我们需要加强测试。
二次探测定理
对于0 < x < p，若p是素数，则方程x^2 \equiv 1 \pmod{p}的解为x_1 = 1, x_2 = p - 1。
证明：原方程可以化为 (x+1)(x-1) \equiv 0 \pmod{p}，即 (x+1)(x-1) 可以整除p。因此x = \pm 1。
因此对于p-1我们考虑拆成d \cdot 2^r。先计算出a^d的值，然后每次平方，检查结果是否为\pm 1，若是，则检查上一次结果是否为\pm 1。如果二次检测未通过，则非素数。这种测试就是Miller-Rabin素性测试了。我们可以知道，这样做的复杂度是O(T \log n)的，T代表测试次数。
但是仍有一些特殊合数能够通过以a为底的MR素性测试，因此我们要不断更换a多次测试。一般来说，我们可以随机几个a出来用，但是也可以预设一些a，比如2 3 7 61 24251这一组能够正确检测大多数数字。可以近似认为选择k个a测试的错误率为0.25^k。
代码
本代码可以通过【P3383】【模板】线性筛素数 – 洛谷一题的较小数据。
// Code by KSkun, 2018/4
#include <cstdio>
#include <cstring>
#include <ctime>
#include <cstdlib>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF
: *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
inline LL fpow(LL n, LL k, LL p) {
LL res = 1;
for(; k; n = n * n % p, k >>= 1) {
if(k & 1) res = res * n % p;
}
return res;
}
inline bool miller_rabin(LL x) {
if(x == 2) return true;
LL t = x - 1, cnt2 = 0;
while(!(t & 1)) {
t >>= 1; cnt2++;
}
for(int i = 0; i < 5; i++) {
LL a = rand() % (x - 1) + 1, now, lst;
now = lst = fpow(a, t, x);
for(int j = 1; j <= cnt2; j++) {
now = lst * lst % x;
if(now == 1 && lst != 1 && lst != x - 1) return false; // 二次检测
lst = now;
}
if(now != 1) return false;
}
return true;
}
LL n, m, t;
int main() {
srand(time(NULL));
n = readint(); m = readint();
while(m--) {
t = readint();
puts(miller_rabin(t) ? "Yes" : "No");
}
return 0;
}
Pollard’s rho算法（Pollard’s rho algorithm）
朴素的分解质因数算法是枚举每个可能的因数，并且把它们全部除出来，这样做的复杂度是O(\sqrt{n} \log n)的。有没有更好的方法？
依赖随机化怎么样？
直接随机到一个因子？概率很小。
那……如果两个因子分别为c和d且c>d，我们让a = \frac{c+d}{2}, b = \frac{c-d}{2}，则a^2 - b^2 = n，枚举a和b，并检查a^2 - n是否是一个完全平方数。概率还是很小啊。
如果我们随机两个数，作差看看是不是呢？emmm没什么差别。
那么如果把作差与n取gcd呢？虽然直接枚举出因子的概率很小，但是因子的倍数是很多的吧。这种想法是能接受的。
于是我们可以随机两个数，作差，如果与n取gcd得到了一个因子则返回它。我们有一个很棒的伪随机函数f(x) = (x^2 + k) \bmod n，其中k是常量，在一个流程中一般不变化，但是取值没有特殊规定。然而，生成的随机数具有循环。
那么我们可以存下哪些数访问过了来避免进入死循环，但是可能会太大存不下呀！于是我们可以利用Floyd判环算法，即设定一个一倍速指针（f(x)），一个两倍速指针（f(f(x))），如果进入环中，两倍速指针一定会在某一时刻追上一倍速指针，此时环肯定已经完整地走完了一圈。（然而底下的代码似乎采用的是算导的方法，即记录下2的幂次结果，判断是否相等）如果遇到环，返回一个寻找失败的信息即可。不要气馁，重新设定常量k再跑一次吧。
其实PR算法的流程就是利用随机化特性找到n的一个因子，然后分治解决。那么分治的终点应该是该数为质数已经不可约了，我们把Miller-Rabin套进去，就可以判断当前数字是否是质数了。
那么总结一下这个算法的流程。
- 如果这个数字已经是质数了，加入这次分解的质因数表中并结束分治。
- 如果非质数，利用PR找到一个因子。
- 将n分为p和n/p两部分，分治解决。
代码
【P1075】质因数分解 – 洛谷
// Code by KSkun, 2018/4
#include <cstdio>
#include <cstring>
#include <ctime>
#include <cstdlib>
#include <vector>
#include <algorithm>
typedef long long LL;
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF
: *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
inline LL fpow(LL n, LL k, LL p) {
LL res = 1;
for(; k; n = n * n % p, k >>= 1) {
if(k & 1) res = res * n % p;
}
return res;
}
LL T, n;
// Miller-Rabin素性检测
inline bool miller_rabin(LL x) {
if(x == 2) return true;
LL t = x - 1, cnt2 = 0;
while(!(t & 1)) {
t >>= 1; cnt2++;
}
for(int i = 0; i < 10; i++) {
LL a = rand() % (x - 1) + 1, now, lst;
now = lst = fpow(a, t, x);
for(int j = 1; j <= cnt2; j++) {
now = lst * lst % x;
if(now == 1 && lst != 1 && lst != x - 1) return false;
lst = now;
}
if(now != 1) return false;
}
return true;
}
std::vector<LL> divi;
inline LL gcd(LL a, LL b) {
if(!b) return a;
return gcd(b, a % b);
}
inline LL randint(LL lim) {
return (1ll * rand() * RAND_MAX + rand()) % lim + 1;
}
inline LL pollard_rho(LL x, LL k) {
LL a = randint(x - 1), b = a, t1 = 1, t2 = 2;
for(;;) {
t1++;
a = (a * a % x + k) % x;
LL g = gcd(abs(b - a), x);
if(g > 1 && g < x) return g;
if(b == a) break;
if(t1 == t2) {
b = a; // 记录2的幂次生成的结果
t2 <<= 1;
}
}
return x;
}
inline void caldivisor(LL x, LL k) {
if(x == 1) return;
if(miller_rabin(x)) {
divi.push_back(x);
return;
}
LL p = x;
while(p >= x) p = pollard_rho(x, k--);
caldivisor(p, k);
caldivisor(x / p, k);
}
int main() {
srand(time(NULL));
n = readint();
caldivisor(n, 200);
printf("%lld\n", *std::max_element(divi.begin(), divi.end()));
return 0;
}
1811 — Prime Test
吐槽：
- 被卡了std::abs的常
- 数据太大，long long也会爆，得换成慢速乘
- ctime头在POJ上会RE
// Code by KSkun, 2018/4
#include <cstdio>
#include <cstring>
//#include <ctime>
#include <cstdlib>
typedef long long LL;
inline LL min(LL x, LL y) {
return x < y ? x : y;
}
inline LL abs(LL x) {
return x < 0 ? -x : x;
}
inline char fgc() {
static char buf[100000], *p1 = buf, *p2 = buf;
return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF
: *p1++;
}
inline LL readint() {
register LL res = 0, neg = 1;
char c = fgc();
while(c < '0' || c > '9') {
if(c == '-') neg = -1;
c = fgc();
}
while(c >= '0' && c <= '9') {
res = res * 10 + c - '0';
c = fgc();
}
return res * neg;
}
inline LL mul(LL a, LL b, LL p) {
a %= p; b %= p;
LL res = 0;
while(b) {
if(b & 1) {
res += a; res %= p;
}
a <<= 1;
if(a >= p) a %= p;
b >>= 1;
}
return res;
}
inline LL fpow(LL n, LL k, LL p) {
LL res = 1;
for(; k; n = mul(n, n, p), k >>= 1) {
if(k & 1) res = mul(res, n, p);
}
return res;
}
LL T, n;
inline bool miller_rabin(LL x) {
if(x == 2) return true;
LL t = x - 1, cnt2 = 0;
while(!(t & 1)) {
t >>= 1; cnt2++;
}
for(int i = 0; i < 20; i++) {
LL a = rand() % (x - 1) + 1, now, lst;
now = lst = fpow(a, t, x);
for(int j = 1; j <= cnt2; j++) {
now = mul(lst, lst, x);
if(now == 1 && lst != 1 && lst != x - 1) return false;
lst = now;
}
if(now != 1) return false;
}
return true;
}
LL mn;
inline LL gcd(LL a, LL b) {
if(a == 0) return 1;
while(b) {
LL t = a % b; a = b; b = t;
}
return a;
}
inline LL pollard_rho(LL x, LL k) {
LL a = rand() % x, b = a, t1 = 1, t2 = 2;
for(;;) {
t1++;
a = (mul(a, a, x) + k) % x;
LL g = gcd(abs(b - a), x);
if(g > 1 && g < x) return g;
if(b == a) return x;
if(t1 == t2) {
b = a;
t2 <<= 1;
}
}
}
inline void caldivisor(LL x) {
if(x == 1) return;
if(miller_rabin(x)) {
mn = min(mn, x);
return;
}
LL p = x;
while(p >= x) p = pollard_rho(x, rand() % (x - 1) + 1);
caldivisor(p);
caldivisor(x / p);
}
int main() {
//srand(time(NULL));
T = readint();
while(T--) {
n = readint();
if(miller_rabin(n)) {
puts("Prime");
continue;
}
mn = 1e15;
caldivisor(n);
printf("%lld\n", mn);
}
return 0;
}
