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

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 时在函数顶部返回 True 来轻松安排= 3，就像你对 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 em>）。所以我们可以立即对其进行测试并为自己节省一次循环迭代：

   # 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 之前，最好先测试可被小素数整除。例如，在Rabin's 1977 paper 他说：

   在实现算法时，我们结合了一些省力的步骤。首先，我们测试任何素数 p N 的可分性，其中 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

Answer 2

除了 Omri Barel 所说的之外，您的 for 循环也存在问题。如果您找到通过测试的a，您将返回true。但是，所有a 都必须通过n 的测试才能成为可能的素数。

Answer 3

我想知道这段代码：

# 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 4

这是我的版本：

# miller-rabin pseudoprimality checker

from random import randrange

def isStrongPseudoprime(n, a):
    d, s = n-1, 0
    while d % 2 == 0:
        d, s = d/2, s+1
    t = pow(a, d, n)
    if t == 1:
        return True
    while s > 0:
        if t == n - 1:
            return True
        t, s = pow(t, 2, n), s - 1
    return False

def isPrime(n, k):
    if n % 2 == 0:
        return n == 2
    for i in range(1, k):
        a = randrange(2, n)
        if not isStrongPseudoprime(n, a):
            return False
    return True

如果你想了解更多关于素数编程的知识，我在我的博客上谦虚地推荐this essay。

Answer 5

您还应该看看Wikipedia，其中已知的“随机”序列给出了给定素数的保证答案。

• 如果n
• 如果n
• 如果 n
• 如果 n
• 如果 n
• 如果 n