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套进去，就可以判断当前数字是否是质数了。
那么总结一下这个算法的流程。

1. 如果这个数字已经是质数了，加入这次分解的质因数表中并结束分治。
2. 如果非质数，利用PR找到一个因子。
3. 将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

吐槽：

1. 被卡了std::abs的常
2. 数据太大，long long也会爆，得换成慢速乘
3. 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;
}

参考资料

• 米勒-拉宾素性检验 – 维基百科，自由的百科全书
• Pollard’s rho algorithm – Wikipedia
• 素数判定的一些讨论（Miller-Rabin算法） – CSDN博客
• poj1811(pollard-rho求一个大数(64位)的最小因子) – CSDN博客
• 大数质因解：浅谈Miller-Rabin和Pollard-Rho算法 – Doggu – 博客园
• 《算法导论（第3版）》，机械工业出版社，9787111407010