Miller–Rabin 质数测试的过程及误差出错概率证明

对于 miller-rabin 的过程,误差和时间复杂度的详细说明。

Miller–Rabin primality test

前言

其实 miller-rabin\text{miller-rabin} 测试我在以前学竞赛的时候就有接触过,而回顾我之前的笔记,发现实际上,有很多地方都未完全弄懂,仅是照搬结论,甚至有错误的部分。因而,这一次,我希望能更深入透彻的了解这个算法。

情景

给一个数 pp,要求判断 pp 是否为质数。这是一个非常常见,也应用十分广泛的问题。而对其的解决方法也值得思考。

首先,一个几乎显然的暴力算法:枚举 2∼p−12\sim p-1 中的所有数判断是否能被 pp 整除。如果不存在能被整除的数,则 pp 是质数。而我们也能稍微优化一下,假设 a∣pa|p,则 pa∣p\frac{p}{a}|p。也就是说,我们枚举的范围可以缩小到 p\sqrt{p}。

当然,即便优化之后,面对很大的质数,这样的时间消耗也是难以接受的。而接下来要介绍的,则是一个比较高效的质数测试算法:miller-rabin\text{miller-rabin}。

费马小定理

在进入 miller-rabin\text{miller-rabin} 之前,首先来介绍一下著名的费马小定理:

若pp为质数,则对于a∈[1,p−1]a\in[1,p-1]有ap≡a( mod p)a^{p}\equiv a(\bmod p)

这里给出我用数学归纳法尝试的一个证明:

显然 a=1a=1 时结论成立。

若 a=na=n 时结论成立,

当 a=n+1a=n+1 时,有:

(n+1)p(n+1)^p=∑i=0pCpinp−i=\sum_{i=0}^{p}C_{p}^{i}n^{p-i}。

那么除了 npn^{p} 和 11 这两项外,

其它的都有一个系数Cpi=p!i!(p−i)!,i∈[1,p−1]C_{p}^{i}=\frac{p!}{i!(p-i)!},i\in[1,p-1]。

由于 pp 为质数,所以 p!p! 中的 pp 并不会在分母中被除掉,所以这些系数都能被pp整除。

而 np≡n( mod p)n^p\equiv n(\bmod p),

所以 (n+1)p≡n+1( mod p)(n+1)^{p}\equiv n+1(\bmod p),结论成立。

因而,当存在 aa 使得 ap−1≡1( mod p)a^{p-1} \not \equiv 1(\bmod p) 时,则 pp 必定不是质数。

至此,可以得到判断 pp 不为质数的一个条件。

可能很自然会想到,如果费马小定理的逆命题成立,我们便能很方便的做出判断。但事实上,有一些合数也能巧妙的满足这样的性质。这类数被称为卡迈克尔数。

这类数的性质非常接近质数,并且证明有无限多个。因此,他们的干扰使得仅靠费马小定理的判断的正确性无法保证。

实际上,我对这类数的条件及性质,以及相关证明也比较感兴趣。但由于本文着重的是 miller-rabin\text{miller-rabin} 相关的方面,便不具体阐述。

miller-rabin\text{miller-rabin}

二次探测定理

费马小定理的局限性使得其会产生失误。那么,可能会自然的想到,能不能加强我们判否的条件,进而使得没有“漏网之鱼”,或者说使其更少呢?

答案是有的,那就是二次探测定理:

对于一个质数 pp, 若 a2≡1( mod p)a^2\equiv1(\bmod p), 则 a≡1( mod p)a\equiv 1 (\bmod p) 或 a≡p−1( mod p)a\equiv p-1(\bmod p)

对其的证明比较简单:

若 a2≡1( mod p)a^2\equiv 1(\bmod p), 则 a2−1≡0( mod p)a^2-1\equiv 0(\bmod p)

即 p∣(a+1)(a−1)p|(a+1)(a-1).

所以只有当 a≡1( mod p)a\equiv 1 (\bmod p) 或 a≡p−1( mod p)a\equiv p-1(\bmod p) 时上式才成立。

所以它也能在 pp 非质数的时候判否,但同样也会有失误而放过合数。

而令我感到有趣的便是,它是如何与费马小定理结合起来的,以及,它能多大程度的加强我们的测试的正确性(或是减小出错的可能性)。

首先,仅考虑 p>2p>2 的情形,则 pp 首先是奇数(否则显然),那么 p−1p-1 就是偶数。

这意味着 ap−1=(ap−12)2a^{p-1}=(a^{\frac{p-1}{2}})^2。

假设 pp 通过了费马小定理的判断,那么 ap−1≡1( mod p)a^{p-1}\equiv1(\bmod p)

也就是说,我们可以对 ap−12a^{\frac{p-1}{2}} 进行二次探测定理的判断。

而且,若其亦为平方数,且 ≡1( mod p)\equiv 1 (\bmod p),则我们又可以对其开方,再继续判断,直至其指数为奇数或 ≡1( mod p)\not\equiv 1(\bmod p)。

在这其中,只要出现不符合二次探测定理的情况,便能立即判否,否则到最后,便有较大概率判断 pp 是质数。

从直观上感觉,我们似乎加上了层层条件,将我们的“网”变得更密了,然而,目前的效果究竟如何?还有没有“鱼”能够漏掉?我们似乎依旧不能确保。接下来,我们就将分析此算法失误的可能性。

失误概率分析

有趣的是,我浏览了目前网络上许多对于 miller-rabin\text{miller-rabin} 的介绍,却并未看到过对其正确性和失误概率的完整证明,包括 wikipedia。大多仅是给出结论:一次测试失误的概率不超过 14\frac{1}{4}。我觉得这可能也是只知其然而不知其所以然的体现,毕竟应用起来确实不需要知道它为什么是对的。

但要深入了解这个算法,则不能避开这个问题。因此,我希望能在这里给出其证明。最终,我在算法导论中找到了对此的证明。然而,这里也仅能证明出错概率不超过 12\frac{1}{2},而更强的结论,需要对需要判断的数组 pp 的选择有所要求,但遗憾的是,我并未找到更加深入探讨这个问题的参考资料。

因此,接下来我将给出关于出错率不超过 12\frac{1}{2} 的证明。


对于数 pp,我们把那些能够在测试中确定 pp 不是质数,即 ap−1≡1( mod p)a^{p-1}\not\equiv 1(\bmod p) 或 a2td≡1( mod p)a^{2^t d}\not\equiv1(\bmod p),的数 aa 称为证据数,反之,则称为非证据数。

则接下来将证明:对于随机选取的数 pp,其非证据数的个数不超过 (p−1)2\frac{(p-1)}{2}。


记 Zp∗\mathbb{Z^*_p} 为模 pp 乘法群(Zp∗,×p)(\mathbb{Z}^*_p,\times_p)的简写。

首先,对于任意非证据数 aa,有 a∈Zp∗a\in \mathbb{Z}_p^*,这是因为

aa 是非证据数

$\Rightarrow a^{p-1}\equiv1(\bmod p) $

⇒\Rightarrow ax≡1( mod p)ax\equiv1(\bmod p) 有解 x=ap−2x=a^{p-2}

⇒gcd⁡(a,p)=1\Rightarrow\gcd(a,p)=1

所以 a∈Zp∗a\in\mathbb{Z}^*_p


取 B={b∈Zp∗∣bp−1≡1( mod p)}B=\{b\in\mathbb{Z}^*_p|b^{p-1}\equiv1(\bmod p)\}

则 BB 在  mod p\bmod p 的乘法下封闭,且对于 ∀\forall 非证据 aa ,有 a∈Ba\in B。那么我们接下来的思路就是限制 BB 的大小。

由群论的拉格朗日定理,

如果HH是一个有限阶的群,而GG是HH的一个子群,那么GG的阶就能整除HH的阶。

如果我们能说明 BB 是 Zp∗\mathbb{Z}^*_p 的一个真子群,则 BB 的大小不超过 p−12\frac{p-1}{2}。接下来,我们将试图证明这一点。而讨论的思路大致如下图所示:

graph LR
p --> A(不是 Carmichael 数)
A --> B(B是真子群)
p --> C(是 Carmichael 数)
C --> D(p 是素数幂)
D --> F(不可能)
C --> E(p 不是素数幂) 
E --> G(B 是真子群)

如果 pp 不是 CarmichaelCarmichael 数,由定义,∃x\exist x 使得 xp−1≡1( mod p)x^{p-1}\not\equiv1 (\bmod p)。也就是说,x∈Zp∗x\in\mathbb{Z}^*_p 且x∈Bx\not\in B。可得 BB 是 Zp∗\mathbb{Z}^*_p 真子群。


否则,接下来证明 pp 不可能是素数幂。

假设 p=qe,e>1p=q^e,e>1,则 qq 是奇素数,由原根定理:

设p是奇素数,则对任意 α\alpha,模 pαp^\alpha 的原根存在。

所以存在原根 gg 使得其阶为 ∣Zp∗∣=ϕ(p)=(q−1)qe−1|\mathbb{Z}^*_p|=\phi(p)=(q-1)q^{e-1}。

又有 gp−1≡1( mod p)g^{p-1}\equiv1(\bmod p),

由离散对数定理:

如果 gg 是 Zn∗\mathbb{Z}^*_n 的一个原根,则当且仅当等式 x≡y( mod ϕ(n))x\equiv y(\bmod \phi(n)) 成立时,有等式 gx≡gy( mod n)g^x\equiv g^y(\bmod n) 成立。

于是,令 y=0y=0 即可得到 ϕ(p)∣p−1\phi(p)|p-1 ,即 (q−1)qe−1∣qe−1(q-1)q^{e-1}|q^e-1。

而此时便产生了矛盾,因为 q∣qe−1q|q^{e-1} 但 q∣qe−1q\not|q^e-1。

所以 pp 不可能是素数幂。

此处的证明可能涉及较多数论的定义以及相关定理,而由于其涉及到的知识较多,又因篇幅关系,不便全部展开一一解释。所以对于此处用到,而能够在他处查证的工具,便没有具体阐述。


最后,我们限制在了 pp 是 CarmichaelCarmichael 数,但不是素数幂的情形下。

若 pp 为合数,则 pp 可分解为 n1×n2n_1\times n_2,其都为奇数且互质。

由算法操作过程,假设 p−1=2t×up-1=2^t\times u,uu 为奇数。

则可以得到一个模 pp 下的序列:X=au,a2u,a22u,…,a2tuX=a^u,a^{2u},a^{2^2u},\dots,a^{2^tu}。

其有两种情况:

  • X=1,1,…,1,1,1,1X=1,1,\dots,1,1,1,1
  • X=x,y,…,−1,1,1,1X=x,y,\dots,-1,1,1,1,即算法到出现 −1-1 时停下,而不管前面的数。

并且,第二种情况一定存在,

取 a=p−1a=p-1,则 (p−1)u≡(−1)u≡−1( mod p)(p-1)^u\equiv(-1)^u\equiv-1 (\bmod p)。

因此,存在一个 vv,使得出现 −1-1 的位置(记为 jj)最靠后,也就是说,v2ju≡−1( mod p)v^{2^ju}\equiv-1(\bmod p),且 j=max⁡k2iu≡−1( mod p)ij=\max\limits_{k^{2^iu}\equiv -1(\bmod p)} i。

则由中国剩余定理的推论:

如果 n1,n2,…,nkn_1,n_2,\dots,n_k 两两互质,且 n=n1n2…nkn=n_1n_2\dots n_k,则对于任意整数 a1,a2,…,aka_1,a_2,\dots,a_k,关于未知量 xx 的联立方程组 x≡ai( mod ni),i=1,2,…,kx\equiv a_i(\bmod n_i),i=1,2,\dots,k对模nn有唯一解。

存在 ww 使得

w≡v( mod n1)w≡1( mod n2)w\equiv v (\bmod n_1) \\ w\equiv 1(\bmod n_2)

由于 v2ju≡−1( mod p)v^{2^ju}\equiv -1 (\bmod p),故有

w2ju≡v2ju≡−1( mod n1)⇒w2ju≡1( mod p)w2ju≡1( mod n2)⇒w2ju≡−1( mod p)w^{2^ju}\equiv v^{2^ju} \equiv -1 (\bmod n_1) \Rightarrow w^{2^ju}\not\equiv 1(\bmod p) \\ w^{2^ju}\equiv 1 (\bmod n_2) \Rightarrow w^{2^ju}\not\equiv -1(\bmod p)

所以 w2ju≡±1w^{2^ju} \not\equiv \pm 1,若 w∈Bw\in B,则在上面的序列中若 w2ku≡−1w^{2^ku}\equiv -1则 k>jk>j ,与 vv 是出现 −1-1 的位置最靠后的而相矛盾。故 w∈Bw\not\in B,则可以得出 BB 是 Zp∗\mathbb{Z}^*_p 的真子集。

至此,我们完成了对于 miller-rabin\text{miller-rabin} 的失误概率分析。那么在足够多次测试后,失误的概率就能被限制到足够小。


从这个证明过程中可以看出,对于其失误概率的分析,却用到了许多如原根,群论等更加深入的数学工具,是将其结合,才得到的这个算法的正确性。但实际上,我自己目前的知识水平有限,对于其中涉及到的知识还了解不深,因此可能目前也还是没有完全吃透。但至少这次,我尝试着把与此算法相关的部分研究了一遍。也算是对自己的一点突破。

而我也不禁好奇,这样的一个算法的提出,究竟是先足够熟练,深刻地了解了各个定理的本质,才足以想到这样的方法,还是先提出了这样的方法后,再用各种工具加以论证?目前的我不得为知,希望今后会有更深的理解。

复杂度分析

同样地,在OIwiki和wikipedia中都仅写了结论,而并未给出相关分析。

在一次测试中,对于 ap−1a^{p-1} 及其之后的二次探测至多共有 log⁡(p−1)\log(p-1) 次,

而每次探测时,若采用快速幂计算 aka^k(kk 为该次的指数),则一次计算为 log⁡k\log k,

同时,对于大数的乘法,需要用类似于快速幂的计算方式,此处有一个 log⁡\log,但亦可以用快速傅里叶变换优化。

同时,为保证正确性,总共进行 TT 次测试。

则复杂度为 Tlog⁡3pT\log^3p。

应用场景

一个比较常见的应用是,在RSA中,需要用到足够大的质数,而此时,用质数筛法等可能耗时较多。而由质数分布,我们知道质数的密度并不算小,因此,我们可以枚举一段连续的数字,而判断其中某个数是否为质数,进而得到我们需要的大质数。利用 miller-rabin\text{miller-rabin},能使这个算法的效率较高。

总结

从最开始的暴力算法,到一步步完善,其实有一个核心的思想:我们并非直接判断 pp 是质数,而是选择判断 pp 不是合数,而更进一步地,将所有 pp 是合数的条件判否。这其中思路的转换,可能有“避其锋芒”的作用,才为整个算法开辟了道路。

而之后,miller-rabin\text{miller-rabin} 体现出了概率算法的巧妙姓。将判断一个数的问题转换为在一个小的概率下失误,而在足够多的尝试后将概率压缩至足够小。实际上,这也是概率算法的魅力所在。而概率算法与确定性算法的关系,至今也是人们还在思考的问题。

同时,一个看上去比较普通的问题,却涉及到如此多较复杂的数学工具,也不得不令人感叹数学的神奇,数字的优美。在此背后,是否蕴含着更加深层的奥秘等待人们揭开?是否有一天,我们能够掌握宇宙的真相,似乎值得思考。

以及,在撰写本文时,阅读我自己的笔记,以及网络上相关资料时,我也发现,似乎并没有哪里,将这个算法讲的足够透彻。无论是 wikipdia 还是 OIwiki,或者是各种人们写的 blog,都总是有令人疑惑而无法在其中得到解答的地方。一方面,可能涉及的东西确实过多,但另一方面,也是否体现人们在研究知识时的局限,仅抓住了核心的结论部分,而忽略了一些证明上,或者别的细节?

同样,本文在某些细节上,可能也没有做到完全的清晰,但我还是希望,能尽量把每个部分讲清楚,能让人在看完后不留疑惑,或者是有更多的兴趣研究更深入的问题。

参考资料

https://www.cnblogs.com/AstatineAi/p/Miller-Rabin-and-Pollard-Rho.html

https://oi-wiki.org/math/number-theory/prime/#miller-rabin-%E7%B4%A0%E6%80%A7%E6%B5%8B%E8%AF%95

https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

https://www.cnblogs.com/zhixingr/p/6750174.html

Thomas H.Cormen.算法导论[M],第三版,殷建平等译,北京:机械工业出版社,2013:567-571.

赞赏