我的素数测试代码有什么问题？

Question

我用Python 3实现了维尔算法（Miller-Rabin素性测试算法），这个算法可以在维基百科上找到。

看起来这个算法在大多数数字上都能正常工作，但偶尔会在某些数字上出错。

比如，素数99999999999999997被判断为不是素数。

我逐行实现了这个算法，但我不知道问题出在哪里。有人能帮我吗？

这是我的代码。

测试输入是：

1

99999999999999997

(两行之间没有空行。)

而我期望的输出应该是“YES”，但在我的机器上却得到了“NO”。

import random

def isPrime(n, k = 5):
'''
Primality test using Miller-Rabin method.
n The number to test primality.
k The number of M-R test to perform.
'''
if n == 1:
    return False
if n == 2 or n == 3:
    return True
if n % 2 == 0:
    return False

# Calculate d
nn = n - 1
s = 1
while nn % (2 ** s) == 0:
    s += 1
s -= 1
d = int(nn / (2 ** s))

for i in range(k):
    a = random.randint(2, n - 1)
    x = pow(a,d,n)
    if x == 1 or x == n - 1:
        continue
    flag = True
    for r in range(1, s):
        x = pow(x,2,n)
        if x == 1:
            return False
        if x == n - 1:
            flag = False
            break
    if not flag:
        continue
    return False
return True

count = int(input())
for i in range(count):
    if isPrime(int(input())):
        print('YES')
    else:
        print('NO')

Answer 1

根据我所了解，Miller-Rabin算法是一种概率算法。你不知道这一点吗？还是说你在用一个修改过的、不是概率性的版本？

Answer 2

我想再说一遍我的评论，因为我的测试结果显示你的例子是有效的。我很怀疑你只是输入测试案例时打错了字。也许你可以再仔细看看？这是我运行后的结果：

在 [12]: millerrabin.isPrime(99999999999999997, 5)

输出 [12]: True

补充：我刚刚运行了更新后的版本，这里是控制台的输出：

1
99999999999999997
YES

再次看，这个结果是正确的。

Answer 3

这是我之前写的一个Miller-Rabin算法的实现。它从来没有给我带来过意外的结果——不过这并不意味着它不会出错！它和你贴的那个基本上是一样的，而且它把99999999999999997判断为素数。我测试你的代码时也是这样判断的——所以这和Mikola的看法一致。但请看下面我发现的一个可能的问题，我一开始没法轻易测试……算了，我测试过了，确实是个问题。

说到素数测试，我不是专家，但我花了很多时间去思考和理解Miller-Rabin算法，我很确定你的实现是正确的。

def is_prime_candidate(self, p, iterations=7):  
    if p == 1 or p % 2 == 0: return False       
    elif p < 1: raise ValueError("is_prime_candidate: n must be a positive integer")
    elif p < self.maxsmallprime: return p in self.smallprimes

    odd = p - 1
    count = 0
    while odd % 2 == 0:
        odd //= 2
        count += 1

    for i in range(iterations):
        r = random.randrange(2, p - 2) 
        test = pow(r, odd, p)
        if test == 1 or test == p - 1: continue
        for j in range(count - 1):
            test = pow(test, 2, p)
            if test == 1: return False
            if test == p - 1: break
        else: return False
        print i
    return True

我注意到你代码中有一点看起来不太对劲：

d = int(nn / (2 ** s))

我心里想，为什么用int呢？然后我意识到你可能是在用Python 3。这就意味着你在这里做的是浮点数运算，然后再转换成整数。这让我觉得有点不靠谱。所以我在ideone上测试了一下。结果是！结果是False。于是我把代码改成了使用明确的向下取整运算（d = nn // (2 ** s)）。结果是！这次是True。