Miller-Rabin算法 - 找不到任何错误？

Question

我用了从Rosetta Code上拿来的代码。我改了一些名字，但其实没做什么大的改动。

import random

def is_probable_prime(n, num_trials = 5):
    assert n >= 2
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    s = 0
    d = n-1
    while True:
        quotient, remainder = divmod(d, 2)
        if remainder == 1:
            break
        s += 1
        d = quotient
    assert(2**s * d == n-1)

    def try_composite(a):
        if pow(a, d, n) == 1:
            return False
        for i in range(s):
            if pow(a, 2**i * d, n) == n-1:
                return False
            return True

    for i in range(num_trials):
        a = random.randrange(2, n)
        if try_composite(a):
            return False
    return True

这个代码和一些伪代码很相似。不过，当我测试这个数字时：

123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901

它返回的是False。而其他用Python和Java写的Miller-Rabin算法实现，测试这个数字时返回的是True，表示它可能是个质数。在经过一些测试后，try_composite在只进行2轮后就返回了True！我真的很想知道哪里出错了，我猜可能是缩进不对或者是我不知道的某个特性。

Answer 1

在你的 try_composite 函数中，for 循环应该写成 for i in range(1,s)。不要测试 i 为零的情况。

编辑：另外，你的 try_composite 函数里缺少一个测试。下面是我写的伪代码版本：

def isPrime(n, k=5):
    def isComposite(s, d):
        x = pow(randrange(2,n-1), d, n)
        if x == 1 or x == n-1: return False
        for r in range(1, s):
            x = pow(x, 2, n)
            if x == 1: return True
            if x == n-1: return False
        return True
    if n < 2: return False
    for p in [2, 3, 5, 7, 11, 13, 17]:
        if n % p == 0: return n == p
    s, d = 0, n-1
    while d % 2 == 0: s, d = s+1, d/2
    for i in range(k):
        if isComposite(s, d): return False
    return True

可惜的是，Python 不允许在 break 或 continue 语句上使用标签。这里有一个更好看的函数伪代码：

function isPrime(n, k=5)
    if n < 2 then return False
    for p in [2,3,5,7,11,13,17,19,23,29]
        if n % p == 0 then return n == p
    s, d = 0, n-1
    while d % 2 == 0
        s, d = s+1, d/2
    for i from 0 to k
        x = powerMod(randint(2, n-1), d, n)
        if x == 1 or x == n-1 then next i
        for r from 1 to s
            x = (x * x) % n
            if x == 1 then return False
            if x == n-1 then next i
        return False
    return True

注意到有两个地方控制流会跳到 next i。在 Python 中没有好的方法来写这个。一个选择是使用一个额外的布尔变量，可以设置和测试来决定何时跳过剩下的代码。另一个选择，就是我上面提到的，写一个局部函数来完成这个任务。这种“循环加一半”的写法很方便也很实用；它在 PEP 3136 中被提出来，但被 Guido 拒绝了。