Rabin-Miller强伪质数测试实现无法正常工作 [英] Rabin-Miller Strong Pseudoprime Test Implementation won't work

问题描述

今天开始尝试实施Rabin-Miller强伪质数测试.

Been trying to implement Rabin-Miller Strong Pseudoprime Test today.

已使用 Wolfram Mathworld 作为参考，第3-5行总结了我的代码很多.

Have used Wolfram Mathworld as reference, lines 3-5 sums up my code pretty much.

但是，当我运行程序时，它说(有时)素数(即使是低的，例如5、7、11)也不是素数.我已经看了很长时间的代码，无法弄清楚出了什么问题.

However, when I run the program, it says (sometimes) that primes (even low such as 5, 7, 11) are not primes. I've looked over the code for a very long while and cannot figure out what is wrong.

为了获得帮助，我在这个站点以及许多其他站点上都进行了研究，但是大多数站点都使用了另一个定义(可能是相同的，但是由于我是这种数学的新手，所以看不到相同的明显联系).

For help I've looked at this site aswell as many other sites but most use another definition (probably the same, but since I'm new to this kind of math, I can't see the same obvious connection).

我的代码:

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给出的定义有何不同?

How does this differ from the definition given at Mathworld?

推荐答案

1.对您的代码进行注释

我在下面将要说明的许多观点在其他答案中都有提及，但将它们结合在一起似乎很有用.

1. Comments on your code

A number of the points I'll make below were noted in other answers, but it seems useful to have them all together.

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:

you've got the roles of r and s swapped: you've actually factored n − 1 as 2^sr. If you want to stick to the MathWorld notation, then you'll have to swap r and s in this section of the code:

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

在线

In the line

for i in range(k):

变量i未使用:通常将此类变量命名为_.

the variable i is unused: it's conventional to name such variables _.

您选择1到 n − 1(含1)之间的随机基数:

You pick a random base between 1 and n − 1 inclusive:

a = random.randrange(1, n)

这就是MathWorld文章中所说的内容，但该文章是从数学家的角度写的.实际上，选择1为底是没有用的，因为1 ^s = 1(mod n )，您会浪费大量的时间进行试验.同样，选择基本 n − 1也没用，因为 s 是奇数，因此( n − 1) ^s = -1(mod n ).数学家不必担心浪费的试验，而程序员则需要担心，因此请写:

This is what it says in the MathWorld article, but that article is written from the mathematician's point of view. In fact it is useless to pick the base 1, since 1^s = 1 (mod n) and you'll waste a trial. Similarly, it's useless to pick the base n − 1, since s is odd and so (n − 1)^s = −1 (mod n). Mathematicians don't have to worry about wasted trials, but programmers do, so write instead:

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

( n 至少需要为4，此优化才能起作用，但是当 n 时，我们可以通过在函数顶部返回True来轻松地安排它= 3，就像 n = 2一样.)

(n needs to be at least 4 for this optimization to work, but we can easily arrange that by returning True at the top of the function when n = 3, just as you do for n = 2.)

如其他答复中所述，您误解了MathWorld文章.当它说" n 通过测试"时，表示" n 通过了基本 a 的测试".关于素数的一个明显事实是，它们通过了 all 个碱基的检验.因此，当您发现 a ^s = 1(mod n )时，您应该做的是四舍五入循环并选择要测试的下一个基准.

As noted in other replies, you've misunderstood the MathWorld article. When it says that "n passes the test" it means that "n passes the test for the base a". The distinguishing fact about primes is that they pass the test for all bases. So when you find that a^s = 1 (mod n), what you should do is to go round the loop and pick the next base to test against.

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

这里有一个优化的机会.我们刚刚计算出的值 x 是 a ^{2 0 s} (mod n ).这样我们就可以立即对其进行测试，并为自己节省一次循环迭代:

There's an opportunity for optimization here. The value x that we've just computed is a^{20 s} (mod n). So we could test it immediately and save ourselves one loop iteration:

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

在计算 a ^{2 j s} 的部分(mod n )中的每个数字都是前一个数字(模 n )的平方.当您只可以平方前一个值时，从头开始计算每个值都是很浪费的.因此，您应将此循环编写为:

In the section where you calculate a^{2j s} (mod n) each of these numbers is the square of the previous number (modulo n). It's wasteful to calculate each from scratch when you could just square the previous value. So you should write this loop as:

# 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

在尝试Miller–Rabin之前，先测试小素数的可除性是一个好主意.例如，在拉宾1977年的论文中，他说:

It's a good idea to test for divisibility by small primes before trying Miller–Rabin. For example, in Rabin's 1977 paper he says:

在实施算法时，我们结合了一些省力的步骤.首先，我们用任何素数 p < N ，例如 N = 1000.

In implementing the algorithm we incorporate some laborsaving steps. First we test for divisibility by any prime p < N, where, say N = 1000.

2.修改后的代码

将所有这些放在一起:

2. Revised code

Putting all this together:

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

这篇关于Rabin-Miller强伪质数测试实现无法正常工作的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！