素性检验

背景

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

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

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

冲突

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

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

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

费马素性检验(Fermat primality test)

费马小定理

若p是素数, a与p互质, 则

维基百科上易懂的证明:

（i）若是整数，是质数，且。若不能整除，则不能整除。取整数集为所有小于的正整数集合（构成的完全剩余系，即中不存在两个数同余），是中所有的元素乘以组成的集合。因为中的任何两个元素之差都不能被整除，所以B中的任何两个元素之差也不能被整除。

换句话说，，考虑共个数，將它们分別除以p，除数分別为，則集合{r₁,r₂,r₃,...,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）。因此

即

在这里W=1·2·3·...·(p*1)，且(W, p) = 1，因此将整个公式除以W即得到：

也即

（ii）若整除，则显然有整除，即。

费马伪素数(Fermat pseudoprime)

对于某些合数x,可以找到与其互质的数a且满足, 则称合数x为费马伪素数, a为x的基. 最小的费马伪素数是341, 以2为基,

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

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

卡迈克尔数(Carmichael number)

卡迈克尔数是费马小定理反面严格不成立的特例. 卡迈克尔数是正合数x, 且使任意与x互质的数a都满足.最小的卡迈克尔数是561.

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

为了检验x是否是素数, 我们可以选择若干个与x互质的数a, 判断是否, 若都满足等于1, 则可以说x有可能是素数.

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

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的解为或

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

证明:

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

设[a, x]已经通过了费马检验, 即, 且x不为偶数(特判2因子).

由于x - 1是偶数, 所以可以转化为二次探测的格式:.

我们的目标是检验x是否是素数, 所以根据二次探测定理, 可以判断是否满足或

这时候会有三种结果:

不符合二次探测定理, 此时可以断定x不为素数
. 如果不是奇数, 则模仿前面的过程, 进行二次探测检验();如果是奇数, 停止检验.
, 停止检验

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

验证

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

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

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')