判断素数的困境
前面的欧拉计划已经遇到很多关于素数的问题。虽然有些题目的核心并非判断素数，但对于那些无法高效解题的题目，如何能够高效判断一个数字是否为素数（又称它为素性测试）的问题便一直萦绕在我心中。传统的试除法是一个确定性算法，然后最终它肯定能给出正确的结果，但时间复杂度最坏可以达到\(Θ(\sqrt{n})\)，在数字很大很大的时候，试除法便会显得十分慢。而筛选法也显得不够经济，虽然优化过的筛选法时间复杂度不高，但是空间复杂度是妥妥的\(O(n)\)，而且往往也是没必要生成这么多素数的。
熟悉JDK的朋友肯定知道，JDK的BigInteger包含了素性测试的方法，主要用在RSA等非对称加密的场景。大数分解、素性测试一直就是现代非对称加密体系的核心。生成非对称密钥对往往会涉及一对十分大的素数，那么JDK究竟用了什么方法能够保证在用户可以接受的短时间内进行素性测试呢？
答案就是米勒-拉宾素性测试。
米勒教授原本是以广义黎曼猜想成立为假设、开发出一个确定性素性测试算法；后来拉宾教授觉得猜想不够靠谱，那个算法就被改成这个不依赖猜想假设的版本。由于我数学较差，也没研究过黎曼猜想，所以就不展开了。刚开始我连这个素性测试算法都看不懂，后来再看了一下它的基础——费马小定理，就基本能明白50%了。后来在一外文网站看到简单的分析（写得比维基百科好的那种），发现一图胜千言，自己画出那个决策图后，便完全明白了。
费马小定理
1636年，业余大数学家费马发现了一个定理：如果p是一个素数，那么对于任意整数a，都有\(a^p – a\)是p的倍数，使用对p取模同余的方式表达，就是\(a^p \equiv a \pmod{p}\)，又或者是常见的\(a^{p-1} \equiv 1 \pmod{p}\)。
对于整数a，如果它是素数p的倍数，即\(a = k \cdot p\)，那么\(a^p – a = [(kp)^{p-1} – 1] \cdot k \cdot p\)，显然是p的倍数，这个毋庸置疑。关键是其它情况的一般化证明。
当然，费马发现它的时候还不是定理，因为没有详细的证明。最早证明的应该是莱布尼茨，人们后来在他1683年前的手稿上发现了这个从未正式发布过的证明。正式发表过证明的是欧拉，在他1736年的《一些与素数有关的定理的证明》论文集中出现。实际上费马小定理是欧拉定理的一个特例，当时欧拉在证明欧拉-费马定理，当欧拉-费马定理成立，费马小定理自然便会成立。
但我不想展开去讨论欧拉-费马定理，因为里面涉及到欧拉函数等等的概念，数学太差看不懂。
目前我能够看得懂的、最直观的费马小定理证明，是利用二项式展开的方式进行的。
首先，对于素数p，假设有一个小于p的正整数n，在组合的运算规则中，从p个当中取出n个，有\(C_p^n = \tbinom{p}{n} = \frac{p!}{n!(p-n)!}\)种选择。由于p是素数，n又小于p，所以分子的p肯定无法约去、即\(\frac{p!}{n!(p-n)!} = p \cdot …\)，即它肯定是p的倍数。所以，对于任意小于p的正整数n，\(\tbinom{p}{n} \equiv 0 \pmod{p}, p > n > 0\)。
而对于任意整数b，关于它的二项式有如下展开过程：
\(\begin{aligned}
(b+1)^p &= C_p^p b^p + C_p^{p-1} b^{p-1} + … + C_p^1 b^1 + C_p^0 b^0 \\
&= \tbinom{p}{p} b^p + \tbinom{p}{p-1} b^{p-1} + … \tbinom{p}{1} b^1 + \tbinom{p}{0} b^0 \\
&= b^p + [\tbinom{p}{p-1} b^{p-1} + … \tbinom{p}{1} b] + 1
\end{aligned}\)
可以观察得到，中间的p-1项，其实都是p的倍数，在对p取模的时候结果都为0，所以便有：\((b+1)^p \equiv b^p + 1 \pmod{p}\)。
然后不断利用二项式的展开，便会有：
\(\begin{aligned}
(b+1)^p &\equiv b^p + 1 \pmod{p} \\
&\equiv [(b-1)+1]^p + 1 \pmod{p} \\
&\equiv [(b-1)^p + 1] + 1 \pmod{p} \\
&\equiv (b-1)^p + 2 \pmod{p} \\
&\equiv [(b-2)+1]^p + 2 \pmod{p} \\
&\equiv [(b-2)^p + 1] + 2 \pmod{p} \\
&\equiv (b-2)^p + 3 \pmod{p} \\
&\equiv … \\
&\equiv [(b-b)+1]^p + b \pmod{p} \\
&\equiv [(b-b)^p+1] + b \pmod{p} \\
&\equiv (b-b)^p + (b + 1) \pmod{p} \\
&\equiv b + 1 \pmod{p}
\end{aligned}\)
然后将a=b+1，直接就可以得到费马小定理的结论\(a^p \equiv a \pmod{p}\)。
能得到这个结论，p是素数是最关键的，因为这直接决定了二项式的系数中的p无法被分母约去，相当于直接将二项式中间的p的倍数的项都在取模运算中被化为0了。这个证明过程不需要用群论、不需要明白欧拉函数，直接有初等数学知识便可以理解。
好了，证明有了，那么费马小定理有什么用？直观地看，它可以用来快速运算指数函数对一个素数取模。然而将费马小定理推广到模不一定是素数（即欧拉-费马定理），只需要模n与底数a互素（互质、即GCD(a, n) = 1），还能够有更广泛的使用空间。这让我联想起《欧拉计划问题48》。
既然费马小定理成立，那么是不是说，我们可以利用它进行素数判断呢？由于a是任意整数，那是不是说，对于整数n，如果有足够多的a使得\(a^n \equiv a \pmod{n}\)成立，便可以称n是一个素数呢？
答案是可以又不可以。因为费马小定理本身只是素数性质的一个必要条件，而不是充分条件——即任何素数都必定满足这个性质，但不是所有满足该性质的数都是素数。数学上将那些能满足该性质的合数称为卡迈克尔数（又称伪素数）。但这不妨碍人们利用费马小定理开发各种各样的素性测试算法，而米勒-拉宾素性测试就是其中一款。
米勒-拉宾素性测试
首先要说明的是，这里对米勒-拉宾素性测试的描述和维基百科里的不一样——因为我数学没那么好、没法子一下子掌握那么多数学术语。
假设p是一个素数，如果有任意整数a，满足\(a^2 \equiv 1 \pmod{p}\)，我们可以很快得到\(a^2 – 1 \equiv 0 \pmod{p}\)，即\((a+1)(a-1) \equiv 0 \pmod{p}\)，即合数(a+1)(a-1)能被素数p整除。由于p是素数，所以要么有(a+1)能被p整除、即\(a+1 \equiv 0 \pmod{p}\)，要么有(a-1)能被p整除、即\(a-1 \equiv 0 \pmod{p}\)。化简：即要么a对p取模余-1、\(a \equiv -1 \pmod{p}\)；要么a对p取模余1、\(a \equiv 1 \pmod{p}\)。
这一组运算有什么用？其实我们可以反过来看，不满足上述情况的数，是不是就可以称它是一个合数？举例我们现在有一个合数c和任意整数a，且满足\(a^2 \equiv 1 \pmod{c}\)。但由于c是合数，那么另一个合数(a+1)(a-1)虽然能被c整除、\((a+1)(a-1) \equiv 0 \pmod{c}\)，但是也无法保证『要么(a+1)能被c整除、要么(a-1)能被c整除』，即无法保证\(a \equiv -1 \pmod{c}\)与\(a \equiv 1 \pmod{c}\)都成立。
也就是说，对于我们要测试的、大于2的奇数n来说，如果能找到一个整数a，能满足\(a^2 \equiv 1 \pmod{n}\)，但\(a \equiv -1 \pmod{n}\)与\(a \equiv 1 \pmod{n}\)都不成立，那我们便可以说n不是一个素数、而是一个合数。
由于n是一个奇数，n-1就肯定是一个偶数。如果n是素数，根据费马小定理，肯定有任意整数a满足\(a^{n-1} \equiv 1 \pmod{n}\)。由于n-1是偶数，肯定可以被2整除，所以根据前面的开方取模运算，便会有要么\(a^{\frac{n-1}{2}} \equiv -1 \pmod{n}\)成立、要么\(a^{\frac{n-1}{2}} \equiv 1 \pmod{n}\)成立。然后观察后面的条件，如果\(\frac{n-1}{2}\)本身依然是一个偶数，那么要判断\(a^{\frac{n-1}{2}} \equiv 1 \pmod{n}\)是否成立，我们可以继续将它展开为判断要么\(a^{\frac{n-1}{4}} \equiv -1 \pmod{n}\)成立、要么\(a^{\frac{n-1}{4}} \equiv 1 \pmod{n}\)。如此类推不断展开，如图所示：
所以实际上，米勒-拉宾素性测试就是不断地将费马小定理以二叉树的方式展开判断的过程。这个多轮判断的过程中，如果n是一个合数/伪素数，是很大机率会被发现的。
为了编程表达方便，我们可以设 q = n-1 ，由于q是一个偶数，所以q又可以表达为\(q = d \cdot 2^k\)，其中d为一奇数。所以实际编程上，我们可以先对\(a^d \equiv 1 \pmod{n}\)进行判断，然后对\(a^d \equiv -1 \pmod{n}\)进行判断、对\(a^{d \cdot 2} \equiv -1 \pmod{n}\)进行判断、对\(a^{d \cdot 2^2} \equiv -1 \pmod{n}\)进行判断、…………最后对\(a^{d \cdot 2^{k-1}} = a^{\frac{q}{2}} = a^{\frac{n-1}{2}} \equiv -1 \pmod{n}\)进行判断。这个判断过程与实际展开过程是相反的，因为这样可以利用编程语言的数学库（或者自己实现）的mod_pow操作简化运算过程。
这个算法的时间复杂度是\(O(i \cdot \log^3{n})\)，i是迭代的次数、n是要判断的奇数。由于每次迭代都要找一个随机数a作为辅助，可以将每次判断都视作一次独立随机事件。单次判断中，假设误判（错把合数当素数、假阳性）的机率不会超过\(\frac{1}{2}\)，也就是说，重复i次判断都误判的机率只有\((\frac{1}{2})^i\)，重复次数越多，出错的可能性越低。不少地方都提及这个算法还有可以提升的空间，我就不继续去搞了，我的目的只是弄明白这套算法的工作原理。
代码实现及测试
结果相当不错，10万范围内的数、每次判断迭代10次，多次测试都没发现误判的情况，平均每次判断1.3毫秒内就可以完成，当数值范围超越JS的内置数值范围，还可以用BigInteger的实现。难怪JDK的BigInteger也优先利用米勒-拉宾算法进行判断，后续再用卢卡斯-莱默算法进行补充。
// 利用新的 W3C API 获取安全的随机数
function rand_uint32() {
const uint32 = new Uint32Array(1)
crypto.getRandomValues(uint32)
return uint32[0]
}
function rand_uint_below(n) {
return rand_uint32() % n + 1
}
function rand_uint_below2(n) {
return Math.ceil(Math.random() * n)
}
// base的平方如果超过JS内置数字范围，这个方法就不可以用了，
// 改为用BigInteger也是可以的，它还自带这个方法
function mod_pow(base, power, mask) {
let prod = 1
for (; power > 0; --power)
prod = prod * base % mask
return prod
}
// 将一个偶数分拆成以2为底数的表达式：even = d*(2^k)，d是一奇数
function split_based_2(even) {
let d = even, k = 0
for (; d % 2 == 0; ++k)
d /= 2
return [d, k]
}
// 利用随机数判断一个奇数是否可能是合数
function is_possible_composite(odd) {
const q = odd - 1 // 对odd取模余数为-1，其实就相当于对odd取模余数为odd-1=q
const a = rand_uint_below(q)
const [d, k] = split_based_2(q) // q是一偶数，q = d*(2^k)
let r = mod_pow(a, d, odd) // 余数 r = a^d % odd
if (r == 1) // a^d ≡ 1 (mod odd)，即odd有可能为素数
return false
if (r == q) // a^d ≡ -1 (mod odd)，即odd有可能为素数
return false
for (let i = 1; i < k; ++i) {
r = r*r % odd // 实际上采用 r = mod_pow(r, 2, odd) 也是可以的
if (r == q) // a^(d*(2^i)) ≡ -1 (mod odd)，即odd有可能为素数
return false
}
return true
}
// Miller-Rabin 素性测试
function is_probable_prime(n, iterations = 10) {
if (n == 1)
return false
else if (n == 2)
return true
else if (n % 2 == 0)
return false
else {
const odd = n
for (let i = 0; i < iterations; ++i)
if (is_possible_composite(odd))
return false
return true
}
}
/**************** TEST ****************/
function* primes_below(n) {
const m = Math.floor((n - 2) / 2)
const sieve = new BitSet
for (let i = 1; i < m; ++i) {
let j = i, p = i + j + 2*i*j
for (; p <= m; ++j, p = i + j + 2*i*j)
sieve.set(p, 1)
}
yield 2
for (let i = 1; i <= m; ++i)
if (sieve.get(i) == 0)
yield (2*i + 1)
}
const limit = 100000
const is_prime = new BitSet
for (const prime of primes_below(limit))
is_prime.set(prime, 1)
function test_miller_rabin() {
try {
const start = Date.now()
for (let i = 0; i < limit; ++i) {
const truth = (is_prime.get(i) == 1)
if (is_probable_prime(i) !== truth) {
if (truth)
// 假阴性
throw new Error('False Negative: ' + i + ' is a prime.')
else
// 假阳性
throw new Error('False Positive: ' + i + ' is a composite.')
}
}
alert('Miller-Rabin: All Good! ' + ((Date.now()-start)/limit) + 'ms/test')
}
catch (err) {
alert(err.toString())
}
}
function test_big_int() {
try {
const start = Date.now()
for (let i = 0; i < limit; ++i) {
const truth = (is_prime.get(i) == 1)
if (bigInt(i).isProbablePrime(10) !== truth) {
if (truth)
// 假阴性
throw new Error('False Negative: ' + i + ' is a prime.')
else
// 假阳性
throw new Error('False Positive: ' + i + ' is a composite.')
}
}
alert('BigInteger: All Good! ' + ((Date.now()-start)/limit) + 'ms/test')
}
catch (err) {
alert(err.toString())
}
}
test_miller_rabin()
test_big_int()
不过这套测试也使我发现了NPM上的BigInteger库的isProbablePrime方法并非采用米勒-拉宾算法的，它实际上是直接使用费马素性测试（写在NPM库页面上）。在10万的范围内，每次判断迭代10次，偶尔会出现假阳性误判——也就是说它显然没处理伪素数的问题！但是如果运气好的话，这套实现可以在更快的速度下完成10万次判断——平均每次不到0.01毫秒！！
由 bruce 于 2018-04-13 23:12 更新
