楕円曲線因数分解アルゴリズムの誕生日のパラドックスの続きを、因数分解プログラムのコレクションに追加したいと思います。ブレントは2つの 論文でアルゴリズムを説明しており、モンゴメリーもアルゴリズムを説明しており、私はBosmaとLenstraによる詳細な説明に従ってアルゴリズムを実装しようとしています。以下は、これまでに Python で作成したもので、 ideone.com /vMXSabで実行できます。
# lenstra's algorithm per bosma/lenstra
from random import randint
from fractions import gcd
def primes(n):
b, p, ps = [True] * (n+1), 2, []
for p in xrange(2, n+1):
if b[p]:
ps.append(p)
for i in xrange(p, n+1, p):
b[i] = False
return ps
def bezout(a, b):
if b == 0: return 1, 0, a
q, r = divmod(a, b)
x, y, g = bezout(b, r)
return y, x-q*y, g
def add(p, q, a, b, m):
if p[2] == 0: return q
if q[2] == 0: return p
if p[0] == q[0]:
if (p[1] + q[1]) % m == 0:
return 0, 1, 0 # infinity
n = (3 * p[0] * p[0] + a) % m
d = (2 * p[1]) % m
else:
n = (q[1] - p[1]) % m
d = (q[0] - p[0]) % m
x, y, g = bezout(d, m)
if g > 1: return 0, 0, d # failure
z = (n*x*n*x - p[0] - q[0]) % m
return z, (n * x * (p[0] - z) - p[1]) % m, 1
def mul(k, p, a, b, m):
r = (0,1,0)
while k > 0:
if k % 2 == 1:
r = add(p, r, a, b, m)
if r[2] > 1: return r
k = k // 2
p = add(p, p, a, b, m)
if p[2] > 1: return p
return r
def lenstra1(n, limit):
g = n
while g == n:
q = randint(0, n-1), randint(0, n-1), 1
a = randint(0, n-1)
b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
g = gcd(4*a*a*a + 27*b*b, n)
if g > 1: return 0, g # lucky factor
for p in primes(limit):
pp = p
while pp < limit:
q = mul(p, q, a, b, n)
if q[2] > 1:
return 1, gcd(q[2], n)
pp = p * pp
return False
def parms(b1):
b2 = 10 * b1
er = [(1,31), (2,63), (3,127), (6,255), (12,511), \
(18,511), (24,1023), (30,1023), (60,2047)]
prev = 1,31
for (e, r) in er:
if e*e > b1/1250: break
prev = e, r
e, r = prev
rBar = int(round(b2/r))
u = randint(0, pow(2,30)//(e+2))
v = randint(0, pow(2,30)//(e+2))
uBar = randint(0, pow(2,30)//(e+2))
vBar = randint(0, pow(2,30)//(e+2))
return b2, e, r, rBar, u, v, uBar, vBar
def lenstra2(n, b1):
g = n
while g == n:
q = randint(0, n-1), randint(0, n-1), 1
a = randint(0, n-1)
b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
g = gcd(4*a*a*a + 27*b*b, n)
if g > 1: return 0, g # lucky factor
for p in primes(b1):
pp = p
while pp < b1:
q = mul(p, q, a, b, n)
if q[2] > 1: return 1, gcd(q[2], n)
pp = p * pp
b2, e, r, rBar, u, v, uBar, vBar = parms(b1)
f = [1] * (r+1)
for i in range(1, r):
p = mul(pow(u*i+v,e), q, a, b, n)
if p[2] > 1: return 2, gcd(p[2], n)
f[i] = (f[i-1] * (q[0] - p[0])) % n
d = 1
for j in range(1, rBar):
pBar = mul(pow(uBar*j+vBar,e), q, a, b, n)
if pBar[2] > 1: return 3, gcd(pBar[2], n)
t = 0
for i in range(0, r):
t = (t + p[0] * f[i]) % n
d = (d * t) % n
g = gcd(d, n)
if 1 < g < n: return 4, g
return False
この関数はエラトステネスのふるいの単純なバージョンを実装し、nprimes未満の素数のリストを返します。関数は拡張ユークリッド アルゴリズムを実装し、aの逆数、bの逆数、およびそれらの最大公約数を返します。楕円演算は関数と関数によって与えられます。add は「ポイント」 (0, 0, d ) を返し、非可逆分母を通知し、それを伝播します。因数分解関数での の使用は、呼び出されるたびにそれをチェックする必要があります。関数は、楕円曲線因数分解の単純な 1 段階バージョンであり、適切に動作します。bezoutaddmulmulmulmullenstra1
関数lenstra2とその補助関数parmsは、上記の Bosma/Lenstra の論文に記載されているアルゴリズムを実装しようとする私の試みです。セクション 6.4 と 6.7 の最適化を考慮せずに、セクション 6.1 で説明されているように、最初に基本的なバージョンを機能させようとしています。parmsの計算は正しいと思います。関数は実行されますが、常に を返しますFalse。これは、因数が見つからなかったことを示します。または、アルゴリズムを完了して最終的なgcd計算から戻る前に、楕円演算の早期中断後に戻ります。問題はfの係数の計算と、 fを使用したdの計算にあると思います。
だから私の質問：
- fの係数を正しく計算しましたか?
- dの値を正しく計算しましたか?
- セクション 6.4 と 6.7 の最適化を実装するにはどうすればよいですか? どちらもわかりません。
- セクション 5.1 の Suyama 曲線を Weierstrass 座標を使用して実装するにはどうすればよいですか?
どうもありがとう。