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 表記法に固執したい場合は、コードのこのセクションでrandを交換する必要があります。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
ラインで
for i in range(k):
変数iは使用されていません: そのような変数に名前を付けるのは慣例です_。
1 からn − 1 までのランダムなベースを選択します。
a = random.randrange(1, n)
これは MathWorld の記事で述べられていることですが、その記事は数学者の視点から書かれています。実際、1 s = 1 (mod n ) であるため、基数 1 を選択しても意味がありません。試行を無駄にすることになります。同様に、基数n − 1 を選ぶのは無意味です。なぜなら、sは奇数であり、したがって ( n − 1) s = −1 (mod n ) だからです。数学者は無駄な試行を心配する必要はありませんが、プログラマーは心配するので、代わりに次のように記述します。
a = random.randrange(2, n - 1)
(この最適化が機能するには、 nは少なくとも 4 である必要がありますが、 n = 2 の場合と同様に、 n = 3Trueの場合は関数の先頭に戻ることで簡単に調整できます。)
他の返信で指摘されているように、あなたは MathWorld の記事を誤解しています。" nがテストに合格する" というのは、" nが基数aのテストに合格する" という意味です。素数に関する際立った事実は、素数がすべての基数のテストに合格することです。したがって、a s = 1 (mod n ) であることがわかったら、ループを一周して、テストする次のベースを選択する必要があります。
# a^(s) = 1 (mod n)?
x = pow(a, s, n)
if x == 1:
continue
ここに最適化の機会があります。計算したばかりの値xは2 0 s (mod n ) です。したがって、すぐにテストして、ループの繰り返しを 1 回節約できます。
# a^(s) = ±1 (mod n)?
x = pow(a, s, n)
if x == 1 or x == n - 1:
continue
a 2 j s (mod n )を計算するセクションでは、これらの各数値は前の数値の 2 乗 (モジュロn ) です。前の値を 2 乗するだけで済むのに、それぞれをゼロから計算するのはもったいないです。したがって、このループを次のように記述する必要があります。
# 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 を試す前に、小さな素数で割り切れるかどうかをテストすることをお勧めします。たとえば、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