素性检验

背景

在知乎上冲浪时了解到一种很有趣的数：其所有非空后缀均为素数的数。例如9613, 其后缀[9613, 613, 13, 3]均为素数. 而此类数的最大值为357686312646216567629137.

证明方法非常直接:用计算机枚举. 不妨为此类数随意的取个名叫做后缀素数, 不难发现, 后缀素数P的后缀必然是后缀素数. 因此对于长度为N的后缀素数必然可以由最高位和一个长度为N-1的后缀素数拼接而成, 如9613可以由是9和613拼接.

因此我们可以从一位的后缀素数集[2, 3, 5, 7]开始, 通过枚举最高位 + 判断是否为素数, 得到下一级位数的所有后缀素数. 若N位的后缀素数集为空, 则最大的后缀素数就存在于N-1位的后缀素数中.

冲突

程序的实现可以分为两部分: 判断素数, 循环递推后缀素数.

程序测试时, 等待了几分种仍然没有结果, 这时我才意识到可能程序的时间复杂度太高.我写的判断素数的函数使用的是非常普通的 $O (\sqrt{n})$ 级别的算法, 在我们已经知道357686312646216567629137为答案的情况下, $O (\sqrt{n})$ 的算法得出结果或许需要几个小时.此时显然需要更高效的素数判断算法.

素性检验是一种概率性检验一个数是否是素数的算法, 换句话说, 素性检验只能得出两种结果: 一个数不是素数或一个数大概率是个素数.

费马素性检验(Fermat primality test)

费马小定理

若p是素数, a与p互质, 则 $a^{p - 1} \equiv 1 (m o d p)$

维基百科上易懂的证明:

（i）若 $a$ 是整数， $p$ 是质数，且 $gcd (a, p) = 1$ 。若 $p$ 不能整除 $x - y$ ，则 $p$ 不能整除 $a (x - y)$ 。取整数集 $A$ 为所有小于 $p$ 的正整数集合（ $A$ 构成 $p$ 的完全剩余系，即 $A$ 中不存在两个数同余 $p$ ）， $B$ 是 $A$ 中所有的元素乘以 $a$ 组成的集合。因为 $A$ 中的任何两个元素之差都不能被 $p$ 整除，所以B中的任何两个元素之差也不能被 $p$ 整除。
换句话说， $gcd (a, p) = 1$ ，考虑 $1 \times a, 2 \times a, 3 \times a, . . . . (p - 1) \times a$ 共 $(p - 1)$ 个数，將它们分別除以p，除数分別为 $r_{1}, r_{2}, r_{3}, . . . ., r_{p - 1}$ ，則集合{r~1~,r~2~,r~3~,...,r~p-1~}为集合{1,2,3,...,(p-1)}的重新排列，即1,2,3,....,(p-1)在除数中恰好各出現一次；这是因为对于任兩个相异k*a而言（k=1,2,3....(p-1)），其差不是p的倍数（所以不会有相同除数），且任一个k*a亦不为p的倍数（所以除数不为0）。因此
$1 \cdot 2 \cdot 3 \cdot \dots \cdot (p - 1) \equiv (1 \cdot a) \cdot (2 \cdot a) \cdot \dots \cdot ((p - 1) \cdot a) (\mod p),$
即
$W \equiv W \cdot a^{p - 1} (\mod p),$
在这里W=1·2·3·...·(p*1)，且(W, p) = 1，因此将整个公式除以W即得到：
$a^{p - 1} \equiv 1 (\mod p)$
也即 $a^{p} \equiv a (\mod p)$
（ii）若 $p$ 整除 $a$ ，则显然有 $p$ 整除 $a^{p}$ ，即 $a^{p} \equiv a \equiv 0 (\mod p)$ 。

费马伪素数(Fermat pseudoprime)

对于某些合数x,可以找到与其互质的数a且满足 $a^{x - 1} \equiv 1 (m o d / x)$ , 则称合数x为费马伪素数, a为x的基. 最小的费马伪素数是341, 以2为基, $2^{340} \equiv 1 (m o d 341)$

对于任意大于1的自然数a, 都有无限个以a为基的费马伪素数.

这也说明, 费马小定理的反面是不成立的, 即使对于x能找到a满足同余式, x也不一定是素数.

卡迈克尔数(Carmichael number)

卡迈克尔数是费马小定理反面严格不成立的特例. 卡迈克尔数是正合数x, 且使任意与x互质的数a都满足 $a^{x - 1} \equiv 1 (m o d x)$ .最小的卡迈克尔数是561.

根据费马小定理, 我们知道, 对于数x, 如果一个数a与x互质, 且 $a^{x - 1} ≢ 1 (m o d x)$ , 则说明x不是一个素数, 将这样的a称为x是合数的凭证(witness).反之, a称为x是素数的强伪证(strong liar).

为了检验x是否是素数, 我们可以选择若干个与x互质的数a, 判断是否 $a^{x - 1} \equiv 1 (m o d x)$ , 若都满足等于1, 则可以说x有可能是素数.

$a$ 通常在区间[2, x-1]中进行选取, 此区间中的数显然都与x互质.

python代码实现如下:

python

def fermat(x, test_time = 8):
    # 小于3的数直接特判
    if x <= 2:
        return x == 2
    for i in range(1, test_time + 1):
        a = randint(2, x - 1)
        # quick_pow为快速幂
        if quick_pow(a, x - 1, x) != 1:
            return False
    return True

米勒-拉宾素性检验(Miller–Rabin primality test)

当检验的数据范围增大时, 费马素性检验的正确率就满足不了人们的需求了.

二次探测定理

如果 $p$ 是奇素数, 则 $x^{2} \equiv 1 (m o d p)$ 的解为 $x \equiv 1 (m o d p)$ 或 $x \equiv - 1 (m o d / p)$ .

特别的, 其小于p的解为 $x \equiv 1 (m o d p)$ 或 $x \equiv p - 1 (m o d p)$

特别提出其小于p的解是因为, 在诸多编程语言中, 模的结果不会取负数.

证明:

x^{2} \equiv 1 (m o d p) ⟺ (x + 1) (x - 1) \equiv 0 (m o d p)

前面我们知道, 在[2, x-1]取一数a, 如果 $a^{x - 1} \equiv 1 (m o d x)$ , 并不能说明x一定是素数, Miller-Rabin素性检验正是在这一步做进一步的纵深, 提高正确率.

设[a, x]已经通过了费马检验, 即 $a^{x - 1} \equiv 1 (m o d x)$ , 且x不为偶数(特判2因子).

由于x - 1是偶数, 所以可以转化为二次探测的格式: $a^{x - 1} \equiv 1 (m o d x) ⟺ (a^{\frac{x - 1}{2}})^{2} \equiv 1 (m o d x)$ .

我们的目标是检验x是否是素数, 所以根据二次探测定理, 可以判断是否满足 $a^{\frac{x - 1}{2}} \equiv 1 (m o d x)$ 或 $a^{\frac{x - 1}{2}} \equiv x - 1 (m o d x)$

这时候会有三种结果:

不符合二次探测定理, 此时可以断定x不为素数
$a^{\frac{x - 1}{2}} \equiv 1 (m o d x)$ . 如果 $\frac{x - 1}{2}$ 不是奇数, 则模仿前面的过程, 进行二次探测检验( $a^{\frac{x - 1}{4}}$ );如果是奇数, 停止检验.
$a^{\frac{x - 1}{2}} \equiv x - 1 (m o d x)$ , 停止检验

python程序实现如下:

python

def miller_rabin(x, test_time = 8):
    if x < 3 | (x & 1):
        return x == 2
    for i in range(1, test_time + 1):
        a = randint(2, x - 1)
        # 费马小定理检验
        if quick_pow(a, x - 1, x) != 1:
            return False
        # 二次探测定理检验
        y = x - 1
        while y % 2 == 0:
            y //= 2
            z = quick_pow(a, y, x)
            if z == x - 1:
                break
            elif z != 1:
                return False
    return True

验证

回到最初的问题上, 现在我们已经有时间复杂度足够低( $l o g^{3} n$ )的算法来判断一个数是否是素数.

后缀素数的递推代码如下:

python

start = time.time()
primes = [2, 3, 5, 7]
j = 10
while True:
    new = []
    for x in primes:
        for i in range(1, 10):
            a = i * j + x
            if (miller_rabin(a)):
                new.append(a)
    if len(new) == 0:
        print("answer:", primes)
        break
    j *= 10
    primes = new
end = time.time()
print('time cost', end - start, 's')