Miller Rabin 算法简介

0.1 一些闲话

最近一次更新是在2019年11月12日。以前的文章有不少问题：当我把个人代码交到LOJ上，发现只有60多分。我调了一个晚上，尝试用{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 61, 24251, 2147483647, 998244353}这么一大串数做为基底，而后左改右改，总算过去了。特别感谢 @骗分过样例 的提醒，如今张贴的代码应该是值得信赖的了。

以前个人同窗好像就指出过个人文章的不少问题。好比说我以前写到，Miller Rabin在面对合数$46856248255981$时，若是用{2,3,7,61,24251}作基底，是会判断错误的。但实际上，他说他对着个人代码写了一遍，发现这个数是能够判掉的。

OI中的数学须要细心。相比其余算法方面，数学真的很差调试——一个公式算错了，一个下标写反了，程序就错了。并且复杂的数学代码很难用gdb查错，只能反复本身的公式是否写对，而且在转换成代码的时候是否有差错。编程须要细致和求实精神。当你在写代码，亦或是在写题解时，多问本身一个：有没有问题？是否是哪里写错了？这里为何要这么写？可不能够造数据Hack？尤为是在写博客的时刻，每个OIer都须要作到足够细致——由于这些文章不是写出来好看，让同窗膜拜的，而是真的要帮到网络另外一端，须要帮助的人。

这是我尤其欠缺的。最近个人博文被指出了不少错误，我也意识到本身的不细致会给本身带来这么多麻烦。也但愿我能做为反面教材，警醒后人吧。

1.1 问题引入

判断一个数$n$是否是质数.

当n不太大时,咱们直接用$O(\sqrt{n})$的试除法就能够了.据说还能够利用素数定理,只用$O(ln(x))$就能够了?这种说法我还暂时没试过.

不过,对于$n \geqslant 10^{18}$时,这个算法是确定不可取的.咱们须要一个更加快速的算法.谁不但愿全部的问题被$O(1)$解决呢?

一想到$O(1)$,咱们不妨试一下随机算法.在以后介绍Pollard Rho算法后,也会用到这种思想.

1.2 随机的依据

按理来说,咱们应该是在$[1,\sqrt{n}]$里面随机,再经过某些方式来提升正确率的.可是这样作太慢!咱们要确保每个在$[1,\sqrt{n}]$中的数都不是n的因子,显然用随机试的算法是行不通的.那么,咱们该怎么办?

数论中,咱们有Fermat(也就是费马)小定理:

$ p \in \text{Prime},\quad x^p \equiv p \pmod{p}$ 或 $ x^{p-1} \equiv 1 \pmod{p}$.

这个定理对于任何质数成立;可是反过来不必定.知足Fermat小定理的数不必定是质数.事实上,若是对任意$ b \in \mathbb{N}, b^{p-1} \equiv 1\pmod{p}$,咱们就称p为Carmichael数.运用Fermat小定理，咱们能够断定一个数不是质数，可是不能断定一个数是质数.

不过,这给咱们的随机枚举提供了一个不错的启示:不是随机枚举$n$的因子,而是一个基底$b$.若对于不一样的基底,都有$b^{p-1} \equiv 1 \pmod{p}$,p是质数的可能性就更大了.事实上,若是p知足$b^{p-1} \equiv 1 \pmod{p}$,那么p就叫作"基于b的伪素数".

若是可能的话,这里会张贴费马小定理的证实.

1.3 算法的改进: 二次探测定理

这个定理叫作"二次"是有缘由的.其一,它是费马小定理探测质数的辅助工具,做为第二道屏障;其二,它和一个"二次同余方程"有关.定理以下:

若$x^2 \equiv 1 \pmod{p}$,则 $x \equiv 1 \pmod{p}$ 或 $x \equiv p-1 \pmod{p}$

这个定理比较好证实.这里是证实过程.

对于一个质数p,只要p和二次探测定理不符,那么咱们就能够确定p不是质数.事实上,咱们经常直接用这个定理去断定质数,而不是用费马小定理.

2.1 基本实现

不管是费马小定理仍是二次探测定理,咱们都须要选取若干个合适的基底来进行测试.通常而言,咱们会选取这五个基底:$2,3,7,61,24251$。当须要测试的数据规模在$10^{16}$左右时，它的表现效果仍是能够的。但当数据规模更大些时，你必须考虑将基底扩大一些。在这个基底的基础上，能够考虑选取前$10\sim 15$大的质数，从而解决几乎$100%$的数据。

对于每个基底$b$和一个待判断的数$x$,咱们进行以下测试: 1.用费马小定理断定一下这个数是否是合数. 2.若是$x-1$是偶数,咱们能够对费马小定理作以下变形:

$b^{x-1}\equiv 1 \pmod{x} \implies b^{\frac{x-1}{2}\cdot 2} \equiv 1 \pmod{x}$ 这样就能够用二次探测定理看一下有没有$b^{\frac{x-1}{2}} \equiv 1 \text{或} x-1 \pmod{x}$ 了.若是都不知足,x就必定不是质数.

若是$\frac{x-1}{2}$一样是一个偶数,咱们能够继续将它拆成$\frac{x-1}{4} \cdot 2$,而后再用一次二次探测定理来作.

3.对于二次探测定理$X^k \equiv 1 \pmod{x}$,若是$k$是奇数或者$X \equiv x-1 \pmod{x}$,此时咱们没有办法再拿这个数作文章了.咱们只能暂时认定$x$是质数.

注意到这就是Miller Rabin算法具备随机性的证据之二:对质数的断定不充分.但事实上,用上面这种方法已经能够成功地断定许多质数了,在实际应用中是值得信赖的.

上面这个算法的时间复杂度是$O((\log n)^2)$的：瓶颈在于快速幂，和分解待验证的数；但实际运行时，速度是至关快的。

(以前这里是有代码的，可是好像有点bug。这里就不张贴了。)

2.2 算法的优化

显然,咱们不难看出上述代码有大量须要优化的地方.就好比说,开始的那次费马小定理的断定是彻底没有必要的.咱们能够直接把这一过程留到二次探测定理去执行. 其次,咱们使用二次探测定理时每次都要探测$\log n$次,这个次数是能够稍微有所优化的.

对于前者,咱们直接删去开头的费马小定理判断就能够了.对于后者,网上广为流传这样一种优化技巧:

设$k=x-1=2^p \cdot d$,$d$是一奇数,那么以前二次探测定理$X^2 \equiv 1 \pmod{x}$检测的数$X$分别是以下几个数:$b^{2^p*d},b^{2^{p-1}*d},b^{2^{p-2}*d},\dots,b^{d}$. 若是咱们从后往前检测,若是其中一个数$X$经过了二次探测定理,就直接断定$x$是质数.

我会尝试着证实这个的.往后可能会贴在这里.

这个新算法的流程以下: 1.按上面的方法计算出d和p.记$cur=b^d$.若是$cur\equiv 1 \text{或} x-1 \pmod{x}$就直接断定x是质数. 2.每次把$cur$赋值为$cur^2 \pmod{x}$直到$cur=b^k$.这个过程等价与将$cur$从$b^d$往$b^k$的方向扫描. 3.若是$cur = x-1$则直接断定x为质数.不然转2. 4.若是上述测试都没有断定x为质数,则直接断定x为合数.

bool mr(ll x,ll b)
{
    ll d,p=0;
    d=x-1;
    while((d&1)==0)
        d>>=1,++p;
    int i;
    ll cur=qp(b,d,x);//快速幂
    if(cur==1 || cur==x-1)
        return true;
    for(i=1;i<=p;++i)
    {
        cur=cur*cur%x; //为了防止溢出,这里最好用快速乘或者直接用__int128转化一下
        if(cur==x-1)
            return true;
    }
    if(i>p)
        return false;
}
bool mr(ll x)
{
    if(x==46856248255981 || x<2)
        return false;
    if(x==2 || x==3 || x==7 || x==61 || x==24251)
        return true;
    return mr(x,2)&&mr(x,3)&&mr(x,7)&&mr(x,61)&&mr(x,24251);
}

算法的具体实现我也是参考了洛谷P4718大佬@Great_Influence 的题解,不太清楚背后的缘由.能够参考一下他的代码.