1

我想将椭圆曲线分解算法的生日悖论延续添加到我的分解程序集合中。Brent 在两篇 论文中描述了该算法，Montgomery 也描述了该算法，我正在尝试根据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

该primes函数实现了一个简单版本的埃拉托色尼筛法，返回一个小于n的素数列表，该bezout函数实现了扩展的欧几里得算法，返回a的倒数、b的倒数以及它们的最大公约数。椭圆算术由addandmul函数给出；add 返回一个“点”(0, 0, d ) 来表示一个不可逆的分母，mul传播它，并且mul在因式分解函数中的使用必须在每次mul调用时检查它。函数lenstra1是椭圆曲线分解的简单单阶段版本，并且可以正常工作。

函数lenstra2及其辅助函数parms是我尝试实现上面引用的 Bosma/Lenstra论文中给出的算法。我首先尝试使基本版本正常工作，如第 6.1 节所述，而不考虑第 6.4 节和第 6.7 节中的优化。我认为里面的计算parms是正确的。函数运行，但总是返回False，表示它没有找到一个因子，或者它在完成算法并从最终gcd计算返回之前在椭圆算术中提前中断后返回。我认为问题在于f的系数的计算，以及使用f来计算d。

所以我的问题：

  1. 我是否正确计算了f的系数？
  2. 我是否正确计算了d的值？
  3. 如何实现 6.4 和 6.7 节的优化？我不明白他们中的任何一个。
  4. 如何使用 Weierstrass 坐标实现第 5.1 节的 Suyama 曲线？

非常感谢。

4

0 回答 0