题解 P4718/论 Miller-Rabin 算法的确定性化

· · 题解

前言

在刚开始学习 Miller-Rabin 的时候，我对算法本身的随机性感到相当程度的不满。

然后，百度百科上的一段话便吸引了我的注意。

卡内基梅隆大学的计算机系教授 Gary Lee Miller 首先提出了基于广义黎曼猜想的确定性算法，
由于广义黎曼猜想并没有被证明，其后由以色列耶路撒冷希伯来大学的 Michael O. Rabin 教授作出修改，
提出了不依赖于该假设的随机化算法。

于是便有了几个问题：

相信不少初学者都有与我一样的问题，
这篇博文将介绍 Miller Rabin 算法及其确定性化，解决这几个疑问。

前置芝士：乘法逆元

费马素性检验

素性检验最暴力的方法是 \Theta(\sqrt{n}) 的，为了应付 2^{64} 以内的数据，我们需要一个更迅速的方法来判定素数。

具体的，我们有费马小定理：

若 p 为素数，对于所有的 1\le a<p，有 a^{p-1}\equiv 1\pmod{p} 。

考虑它的逆否命题：

若存在 1\le a<p，使得 a^{p-1}\not\equiv 1\pmod{p}，则 p 一定不是素数。

那我们可以在 [1,p-1) 内随便选取几个数，快速幂来检验 p 是否是素数。

这一方法称为 费马素性检验 。

二次探测定理与 Miller Rabin 素性检验

遗憾的是，费马素性检验存在一个问题：
有的合数 p 也满足 a^{p-1}\equiv 1\pmod{p} \quad(1\le a<p)，
这样的数被称为卡迈克尔数（Carmichael Number），又称为费马伪素数。
例如，561=3\times 11\times 17 就是一个卡迈克尔数。
在 10^9 以内，这样的数有 646 个，显然 2^{64} 以内不能通过打表来满足素性检测的要求。

这迫使我们去优化费马素性检验。具体地，我们有二次探测定理：

若 p 为奇素数，则 x^2\equiv 1\pmod{p} 的解为 x\equiv\pm 1\pmod{p}。

若 p 是奇数，显然 p-1 是偶数，可以这样优化算法的正确性：

typedef unsigned long long ull;
typedef unsigned int word;
bool check(const word a,const ull p){
    ull d=p-1,get=pow(a,d,p);
    if(get!=1) return 1;//特判 d=p-1 的情况
    while((d&1)^1)
        if(d>>=1,(get=pow(a,d,p))==p-1) return 0;//先 d/=2，再计算快速幂
        else if(get!=1) return 1;
    return 0;
}

初学者写 Miller-Rabin 时可能会犯的错误：

若随机选取 k 个底数，Miller-Rabin 的错误率为 4^{-k}。
也就是说，选取 20\sim 40 个底数，正确率是可以接受的。

固定底数与确定性检验

Miller-Rabin 的时间复杂度为 \Theta(k\log^2 p)，实际是用于获取高达 2^{1024} 的 “工业” 素数的。
在 OI 的 2^{64} 范围内，我们可以将其确定性化。

如果我们假设广义 Riemann 猜想（GRH）成立，那么证明 n 是素数就只需要验证
“ n 可以通过以 2,3,4,\cdots \lfloor 2(\log n)^2\rfloor 为底的 Rabin-Miller 测试 ”了。
这个算法，如果能严格证明的话，运行时间是 O((\log n)^5)。
——摘自知乎回答如何编程判断一个数是否是质数？

即使 GRH 成立，改造后算法的运行时间已经不可接受了。 但是，就算 GRH 不成立，通过选取几个固定的底数， 我们依然能在一定的范围内实现确定性判素。 - 对于 $2^{32}$ 以内的判素，选取 $2,7,61$ 三个底数即可。 - 对于小于 $2^{64}$ 以内的判素，选取 $2,325,9375,28178,450775,9780504,1795265022$ 七个底数即可。 - 如果是考场上，选取 $2,3,5,7,11,13,17,19,23,29,31,37

（也就是前十二个质数）作为底数即可，它适用于 2^{78} 以内的判素。
对于选取前 n 个质数作为底数的适用范围，详见 OEIS

由此我们得到了最终的固定判素的代码（需要 C++ 11）：

typedef unsigned long long ull;
typedef unsigned char byte;
const byte test[]={2,3,5,7,11,13,17,19,23,29,31,37};
bool miller_rabin(const ull p){
    if(p>40){
        for(word a:test)
            if(check(a,p)) return 0;
        return 1;
    }
    for(word a:test)
        if(p==a) return 1;
    return 0;
}

参考资料

希望借我这篇博文，确定性判素的 miller-rabin 能在 OI 中普及。
望大家学习时不要人云亦云、浅尝辄止。