<h3>一。对代码的评论</h3>
<p>下面我要说的几点在其他的答案中也提到了,但是把它们放在一起似乎很有用。</p>
<ol>
<li><p>在本节中</p>
<pre class="lang-py prettyprint-override"><code>s = 0
r = n - 1
# factor n - 1 as 2^(r)*s
while r % 2 == 0:
s = s + 1
r = r // 2 # floor
</code></pre>
<p>你已经交换了<em>r</em>和<em>s</em>的角色:你实际上已经将<em>n</em>-1分解为2<sup><em>s</em></sup><em>r</em>。如果要坚持使用MathWorld表示法,则必须在代码的这一部分中交换<code>r</code>和<code>s</code>:</p>
<pre class="lang-py prettyprint-override"><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
</code></pre></li>
<li><p>排队</p>
<pre><code>for i in range(k):
</code></pre>
<p>变量<code>i</code>未使用:通常将此类变量命名为<code>_</code>。</p></li>
<li><p>在1和<em>n</em>-1之间选择一个随机基数:</p>
<pre><code>a = random.randrange(1, n)
</code></pre>
<p>这是它在MathWorld文章中说的,但那篇文章是从数学家的角度写的。事实上,选择基数1是没用的,因为1<sup><em>s</em></sup>=1(mod<em>n</em>)会浪费一次试验。类似地,选择基<em>n</em>-1也没用,因为<em>s</em>是奇数,所以(<em>n</em>-1)<sup><em>s</em></sup>=-1(mod<em>n</em>)。数学家不必担心试验的浪费,但是程序员需要,所以写下:</p>
<pre><code>a = random.randrange(2, n - 1)
</code></pre>
<p>(<em>n</em>要使此优化工作,至少需要4,但是我们可以在<em>n</em>=3时在函数顶部返回<code>True</code>来轻松地安排,就像您对<em>n</em>=2所做的那样)</p></li>
<li><p>如其他回复所述,您误解了MathWorld文章。当它说“<em>n</em>通过测试”时,它意味着“<em>n</em>通过了基础<em>a</em>的测试”。素数的区别在于它们通过了对所有碱基的检验。因此,当你发现<em>a</em><sup><em>s</em></sup>=1(mod<em>n</em>)时,你应该做的是绕过这个循环并选择下一个要测试的基。</p>
<pre><code># a^(s) = 1 (mod n)?
x = pow(a, s, n)
if x == 1:
continue
</code></pre></li>
<li><p>这里有一个优化的机会。我们刚刚计算的值是<em>a</em><sup>2<sup>0<sup>s<sup>(mod<em>n</em>)。因此,我们可以立即测试它,并为自己节省一个循环迭代:</p>
<pre><code># a^(s) = ±1 (mod n)?
x = pow(a, s, n)
if x == 1 or x == n - 1:
continue
</code></pre></li>
<li><p>在计算<em>a</em><sup>2<sup><em>j</em></sup>s</sup>(mod<em>n</em>)的部分中,每个数字都是前一个数字的平方(模<em>n</em>)。当你能把前一个值平方时,从头计算每一个值是浪费的。所以你应该把这个循环写成:</p>
<pre><code># 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
</code></pre></li>
<li><p>在尝试Miller-Rabin之前,最好先测试一下小素数的可分性。例如,在<a href="http://www.sciencedirect.com/science/article/pii/0022314X80900840" rel="nofollow">Rabin's 1977 paper</a>中,他说:</p>
<blockquote>
<p>In implementing the algorithm we incorporate some laborsaving steps. First we test for divisibility by any prime <em>p</em> < <em>N</em>, where, say <em>N</em> = 1000.</p>
</blockquote></li>
</ol>
<h3>2。修订代码</h3>
<p>把这些放在一起:</p>
<pre class="lang-py prettyprint-override"><code>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
</code></pre>