Frobenius 疑似素数アルゴリズムは、Miller-Rabin 素数性テストよりも実行に 3 倍の時間がかかりますが、解像度は 7 倍になると誰かが私に言いました。したがって、前者を 10 回実行し、後者を 30 回実行すると、両方の実行に同じ時間がかかりますが、前者は約 233% 高い分析能力を提供します。テストを実行する方法を見つけようとしているときに、最後にアルゴリズムが記載された次の論文が発見されました。
以下のアルゴリズムを実装しようとしていますが、プログラムは数値を出力しません。数学表記法やアルゴリズムに詳しい人は、何が起こっているのかを確認してもらえますか?
編集 1:以下のコードには修正が追加されていますが、実装がcompute_wm_wm1ありません。誰かがアルゴリズムの観点から再帰的な定義を説明できますか? 私にとっては「クリック」ではありません。
編集 2:誤ったコードが削除され、compute_wm_wm1関数の実装が以下に追加されました。機能しているように見えますが、実用化するにはさらに最適化が必要になる場合があります。
from random import SystemRandom
from fractions import gcd
random = SystemRandom().randrange
def find_prime_number(bits, test):
number = random((1 << bits - 1) + 1, 1 << bits, 2)
while True:
for _ in range(test):
if not frobenius_pseudoprime(number):
break
else:
return number
number += 2
def frobenius_pseudoprime(integer):
assert integer & 1 and integer >= 3
a, b, d = choose_ab(integer)
w1 = (a ** 2 * extended_gcd(b, integer)[0] - 2) % integer
m = (integer - jacobi_symbol(d, integer)) >> 1
wm, wm1 = compute_wm_wm1(w1, m, integer)
if w1 * wm != 2 * wm1 % integer:
return False
b = pow(b, (integer - 1) >> 1, integer)
return b * wm % integer == 2
def choose_ab(integer):
a, b = random(1, integer), random(1, integer)
d = a ** 2 - 4 * b
while is_square(d) or gcd(2 * d * a * b, integer) != 1:
a, b = random(1, integer), random(1, integer)
d = a ** 2 - 4 * b
return a, b, d
def is_square(integer):
if integer < 0:
return False
if integer < 2:
return True
x = integer >> 1
seen = set([x])
while x * x != integer:
x = (x + integer // x) >> 1
if x in seen:
return False
seen.add(x)
return True
def extended_gcd(n, d):
x1, x2, y1, y2 = 0, 1, 1, 0
while d:
n, (q, d) = d, divmod(n, d)
x1, x2, y1, y2 = x2 - q * x1, x1, y2 - q * y1, y1
return x2, y2
def jacobi_symbol(n, d):
j = 1
while n:
while not n & 1:
n >>= 1
if d & 7 in {3, 5}:
j = -j
n, d = d, n
if n & 3 == 3 == d & 3:
j = -j
n %= d
return j if d == 1 else 0
def compute_wm_wm1(w1, m, n):
a, b = 2, w1
for shift in range(m.bit_length() - 1, -1, -1):
if m >> shift & 1:
a, b = (a * b - w1) % n, (b * b - 2) % n
else:
a, b = (a * a - 2) % n, (a * b - w1) % n
return a, b
print('Probably prime:\n', find_prime_number(300, 10))