Rabin-Miller Strong Pseudoprime Test Implementation将无效

Question

今天一直试图实施Rabin-Miller Strong Pseudoprime Test.

使用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 ^s r.如果你想坚持到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才能使此优化起作用,但我们可以通过True在n = 3 时返回函数的顶部来轻松排列,就像你对n = 2一样.)

4. 正如其他回复中所述,您误解了MathWorld的文章.当它说" n通过测试"时,意味着" n通过测试基础a ".关于素数的显着特征是它们通过了所有基础的测试.所以,当你发现一个^小号 = 1(MOD ñ),你应该做的是去一轮循环,并挑选下一个基地要测试的.

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

5. 这里有优化的机会.值X,我们刚刚计算是一个^{2 0小号}(MOD ñ).所以我们可以立即测试它并保存自己一个循环迭代:

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

6. 在你计算部一个^{2 Ĵ小号}(MOD Ñ),每个这些数字是以前的数(模数的平方Ñ).当你可以将前一个值平方时,从头开始计算它们是浪费的.所以你应该把这个循环写成:

   # 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. 在尝试米勒 - 拉宾之前,测试小素数的可分性是个好主意.例如,在Rabin 1977年的论文中,他说:

   在实现算法时,我们采用了一些省力的步骤.首先,我们通过任何素数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