Rabin-Miller强伪素数测试实现无法工作

Question

今天我在尝试实现拉宾-米勒强伪素数测试。

我参考了Wolfram Mathworld，其中第3到第5行基本上总结了我的代码。

但是，当我运行程序时，有时它会说一些素数（甚至像5、7、11这样的小素数）不是素数。我仔细检查了代码很久，但还是找不到问题出在哪里。

为了寻求帮助，我查看了这个网站以及其他很多网站，但大多数使用的是另一种定义（可能是相同的，但因为我对这种数学还很陌生，所以无法看到明显的联系）。

我的代码：

import random

def RabinMiller(n, k):

    # obviously not prime
    if n < 2 or n % 2 == 0:
        return False

    # special case        
    if n == 2:
        return True

    s = 0
    r = n - 1

    # factor n - 1 as 2^(r)*s
    while r % 2 == 0:
        s = s + 1
        r = r // 2  # floor

    # k = accuracy
    for i in range(k):
        a = random.randrange(1, n)

        # a^(s) mod n = 1?
        if pow(a, s, n) == 1:
            return True

        # a^(2^(j) * s) mod n = -1 mod n?
        for j in range(r):
            if pow(a, 2**j*s, n) == -1 % n:
                return True

    return False

print(RabinMiller(7, 5))

这和Mathworld给出的定义有什么不同呢？

Answer 1

我在想这段代码：

# factor n - 1 as 2^(r)*s
while r % 2 == 0:
    s = s + 1
    r = r // 2  # floor

假设我们有 n = 7。那么 n - 1 = 6。我们可以把 n - 1 表示成 2^1 * 3。在这个情况下，r = 1，s = 3。

但是上面的代码找到了其他的东西。它开始时设定 r = 6，所以 r % 2 == 0。最开始 s = 0，经过一次循环后，s = 1，而 r = 3。但这时 r % 2 != 0，循环就结束了。

最后我们得到 s = 1 和 r = 3，这显然是不对的：2^r * s = 8。

你不应该在循环中更新 s。相反，你应该计算能把数字除以 2 的次数（这就是 r），而经过这些除法运算后的结果就是 s。在 n = 7 的例子中，n - 1 = 6，我们可以除一次（所以 r = 1），除完后得到 3（所以 s = 3）。

Answer 2

除了Omri Barel说的内容，还有一个问题在于你的for循环。你只要找到一个通过测试的就会返回true。但是，要想让n被认为是一个可能的质数，所有的都必须通过测试。

Answer 3

1. 代码中的注释

下面我提到的一些要点在其他回答中也有提到，但把它们集中在一起看会更有帮助。

1. 在这一部分

   s = 0
   r = n - 1
   
   # factor n - 1 as 2^(r)*s
   while r % 2 == 0:
       s = s + 1
       r = r // 2  # floor
   

   你把 r 和 s 的角色搞混了：实际上你应该把 n − 1 分解成 2^sr。如果你想保持 MathWorld 的符号，那么在这段代码中你需要交换 r 和 s：

   # factor n - 1 as 2^(r)*s, where s is odd.
   r, s = 0, n - 1
   while s % 2 == 0:
       r += 1
       s //= 2
   

2. 在这一行

   for i in range(k):
   

   变量 i 没有被使用：通常这种变量会命名为 _。

3. 你在 1 到 n − 1 之间随机选择一个基数：

   a = random.randrange(1, n)
   

   这在 MathWorld 的文章中是这么说的，但那篇文章是从数学家的角度写的。实际上，选择基数 1 是没用的，因为 1^s = 1 (mod n)，这样你就浪费了一次测试。同样，选择基数 n − 1 也是没用的，因为 s 是奇数，所以 (n − 1)^s = −1 (mod n)。数学家不需要担心浪费测试，但程序员需要，所以你应该改成：

   a = random.randrange(2, n - 1)
   

   (n 至少需要为 4 才能让这个优化生效，但我们可以很容易地处理这个问题，当 n = 3 时直接返回 True，就像你对 n = 2 时做的那样。)

4. 正如其他回复中提到的，你误解了 MathWorld 的文章。当它说 "n 通过了测试" 时，它的意思是 "n 对基数 a 通过了测试"。质数的一个重要特征是它们对 所有 基数都通过测试。所以当你发现 a^s = 1 (mod n) 时，你应该继续循环，选择下一个基数进行测试。

   # a^(s) = 1 (mod n)?
   x = pow(a, s, n)
   if x == 1:
       continue
   

5. 这里有一个优化的机会。我们刚计算出的值 x 是 a^{20 s} (mod n)。所以我们可以立即测试它，省去一次循环：

   # a^(s) = ±1 (mod n)?
   x = pow(a, s, n)
   if x == 1 or x == n - 1:
       continue
   

6. 在你计算 a^{2j s} (mod n) 的部分，每个数字都是前一个数字的平方（模 n）。从头开始计算每一个是浪费的，你可以直接平方前一个值。所以你应该把这个循环写成：

   # a^(2^(j) * s) = -1 (mod n)?
   for _ in range(r - 1):
       x = pow(x, 2, n)
       if x == n - 1:
           break
   else:
       return False
   

7. 在尝试 Miller–Rabin 之前，先测试一下能否被小质数整除是个好主意。例如，在 拉宾的 1977 年论文中他说：

   在实现算法时，我们加入了一些省力的步骤。首先，我们测试是否能被任何小于 N 的质数 p 整除，比如说 N = 1000。

2. 修订后的代码

把这些内容整合在一起：

from random import randrange

small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] # etc.

def probably_prime(n, k):
    """Return True if n passes k rounds of the Miller-Rabin primality
    test (and is probably prime). Return False if n is proved to be
    composite.

    """
    if n < 2: return False
    for p in small_primes:
        if n < p * p: return True
        if n % p == 0: return False
    r, s = 0, n - 1
    while s % 2 == 0:
        r += 1
        s //= 2
    for _ in range(k):
        a = randrange(2, n - 1)
        x = pow(a, s, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True