Rabin Miller强伪素数测试实现

Question

今天一直在尝试实现Rabin-Miller强伪素数测试。

使用了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循环也有问题。如果找到一个通过测试的a，则返回true。然而，所有的a都必须通过n测试才能成为可能的素数。

Answer 3

一。对代码的评论

下面我要说的几点在其他的答案中也提到了，但是把它们放在一起似乎很有用。

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（modn）会浪费一次试验。类似地，选择基n-1也没用，因为s是奇数，所以（n-1）^s=-1（modn）。数学家不必担心试验的浪费，但是程序员需要，所以写下：

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

   （n要使此优化工作，至少需要4，但是我们可以在n=3时在函数顶部返回True来轻松地安排，就像您对n=2所做的那样）

4. 如其他回复所述，您误解了MathWorld文章。当它说“n通过测试”时，它意味着“n通过了基础a的测试”。素数的区别在于它们通过了对所有碱基的检验。因此，当你发现a^s=1（modn）时，你应该做的是绕过这个循环并选择下一个要测试的基。

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

5. 这里有一个优化的机会。我们刚刚计算的值是a^{20s（modn）。因此，我们可以立即测试它，并为自己节省一个循环迭代：}

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

6. 在计算a^2js（modn）的部分中，每个数字都是前一个数字的平方（模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之前，最好先测试一下小素数的可分性。例如，在Rabin's 1977 paper中，他说：

   In implementing the algorithm we incorporate some laborsaving steps. First we test for divisibility by any prime p < N, where, say 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