伦斯特拉椭圆曲线分解问题

Question

我正在尝试使用Hendrik Lenstra的椭圆曲线分解方法来分解小的(小于40位)复合整数。

import math 
from fractions import gcd
import random 
def lenstra_elliptic_curve_factor(N):
""" Lenstra's elliptic curve factoring method """
if N <=0:
raise Exception("Integer %s must be possitive " % N) 
# Can't be 1 and can't factor a prime! 
if 1 <= N <= 2 or is_probable_prime(N):
return [N]
# random point in the plain (values less than N)
x0, y0 = random.randrange(1, N), random.randrange(1, N)
factors = list()
bound = int(math.sqrt(N))
for a in xrange(2,N):
# Build curve out of random points
b = y0**2 - x0**3 - a*x0
# Check curve is not singular 
if 4*a**3 - 27*b**2 ==0:
continue
# Initially double point 
s = 3*x0**2 + a
(x,y) = (s**2 - 2*x0, s*((s**2 - 2*x0) - x0) - y0)
# Keep adding points until gcd(x-x0,N) != 1
for k in xrange(2,bound):
for i in xrange(0,math.factorial(k)):
d = gcd(x- x0,N)
if d != 1:
return lenstra_elliptic_curve_factor(int(d)) + lenstra_elliptic_curve_factor(int(N/d))
else:
# Point doubling arithmetic 
s = (y - y0) * modInv(x - x0,N)
x = s**2 - x - x0  
y = - y + s * (s**2 - x - x0 - x)

其中is_probably_prime()是米勒-拉宾检验，试验次数设置为 20。似乎对于某些合数，例如 4，它没有找到非平凡gcd(x-x0)，相反，算法一直通过并且什么也不返回。因此，当算法尝试分解 4 除以的较大数字时，例如 12，return lenstra_elliptic_curve_factor(int(d)) + lenstra_elliptic_curve_factor(int(N/d))向列表中添加一个"NoneType"。例如

for x in xrange(0,3241):
print x, lenstra_elliptic_curve_factor(x) 

我得到

0 [0]
1 [1]
2 [2]
3 [3]
4 None
5 [5]
6 None
7 [7]
8 None
9 [3, 3]
10 [2, 5]
11 [11]
12
Traceback (most recent call last):
File "/AVcrypto/util.py", line 160, in <module>
print x, lenstra_elliptic_curve_factor(x) 
File "/Users/Kevin/gd/proj/curr/honproj/AVcrypto/util.py", line 104, in lenstra_elliptic_curve_factor
return lenstra_elliptic_curve_factor(int(d)) + lenstra_elliptic_curve_factor(int(N/d))
TypeError: can only concatenate list (not "NoneType") to list

我尝试将测试的曲线数量增加到N**10但似乎具有相同的结果。我只是想知道是否有人对这种算法有任何经验，特别是某些数字似乎在很长一段时间内避免了审判过程。

Answer 1

Lenstra 的算法假设被分解的数字是 6 的互质数(即没有 2 或 3 的因数)。尝试因素 4 是行不通的。更现实的测试是考虑13290059因素。

我假设您知道 ECM 对于 40 位数字来说是巨大的矫枉过正。

Answer 2

有趣的问题，谢谢！

决定从头开始实现我自己的ECM方法版本。它有相当小且编程干净的代码，所以希望它对你有用，弄清楚如何改进你的代码。

关于你的代码。一个主要缺点是你没有创建很多曲线，你只使用了一条曲线。经典 ECM 生成许多具有增长边界的曲线。第二个较小的缺点是使用阶乘作为曲线乘数，更快的是使用小素数的小幂乘积，它就足够了，而不是阶乘。

代码中非常重要的缺点是，您确实对自身进行了阶乘时间的点加法。如果您阅读文章椭圆曲线乘法，您会发现对于乘法，执行二进制幂就足够了，就像对常规数字所做的一样，如从右到左的二进制幂中所述。这种幂比你的电流快得多，尤其是当阶乘变大时。现在你做O(Factorial)加法，而你可以只做O(Log2(Factorial))，少得多的操作。

您还说代码的主要问题是它跳过了非常小的因素，例如 4。之所以会这样，是因为 ECM 适用于更大的数字。对于几个微小的数字，如4，它可能有一些特殊的行为。在这种情况下，常规做法是使用试验除法分解来删除所有较小的因素，例如 16 位。对于所有大于 16 位的数字，ECM 肯定会起作用。

我完全从头开始实现了所有(甚至是小的)子算法，甚至没有使用任何库函数 math.gcd()。库外函数仅使用 math.sqrt() 和 random.randrange()。

我在代码中使用以下子算法：试验除法分解、费马概率检验、埃拉托色尼筛(素数生成器)、欧几里得算法(计算最大公约数，GCD)、扩展欧几里得算法(带有 Bezu 系数的 GCD)、模块化乘法逆、从右到左的二进制幂(用于椭圆点乘法)、椭圆曲线算术(点加法和乘法)、Lenstra 椭圆曲线因式分解。

我根据 Wiki 中的 ECM 的算法如下：

1. 检查数字是否小于 2^16，然后通过试除法将其分解。返回结果。

2. 检查数字是否可能是素数，具有高条件，为此我使用费马测试进行 32 次试验。为了克服卡迈克尔数字，你可以使用米勒拉宾测试代替。如果数字是素数，则将其作为唯一的因子返回。

3. 随机生成曲线参数A, X, Y并从曲线方程Y^2 = X^3 + AX + B (mod N)推导出B。检查曲线是否正常，值4 * A ** 3 - 27 * B ** 2应为非零。

4. 通过埃拉托色尼筛生成小素数，素数低于我们的界限。每个素数提高到某个小幂，这个提高的素数将称为K。我确实在它比一些 Bound2 小的时候提高到功率，这在我的情况下Sqrt(Bound)。

5. 计算椭圆点乘法P = k * P，其中 K 取自上一步，P 为 (X， Y)。根据维基计算。

6. 点乘法使用模逆，它根据 Wiki 计算GCD(SomeValue, N)。如果此 GCD 不是 1，则它给出 N 的非 1 因子，因此在本例中 I 通过异常并从 ECM 分解算法返回此因子。

7. 如果所有素数直到 Bound 相乘并且没有给出因子，则使用另一条随机曲线和更大的 Bound 再次重新运行 ECM 分解算法(1.-6.上)。在我的代码中，我通过将 256 添加到旧边界来获取新绑定。

作为以下代码的示例，我在运行后 0.1-5 秒内分解 72 位数字，非常快！多次重新运行代码以查看不同示例数字的分解。如果花费太多时间，则中断程序并重新运行它以随机推出另一个更简单的数字。

96 位(29 位)随机数的控制台输出示例：

N: 41498298392926620982104401176
Curve    0,  bound 2^  8.00
Factor ECM: 8
Factors TrialDiv: [2, 2, 2]
Curve    1,  bound 2^  9.00
Factor ECM: 13
Curve    2,  bound 2^  9.58
Factor ECM: 11
Curve    3,  bound 2^ 10.00
Curve    4,  bound 2^ 10.32
Factor ECM: 45061
Curve    5,  bound 2^ 10.58
Curve    6,  bound 2^ 10.81
Curve    7,  bound 2^ 11.00
Factor ECM: 2338703
41498298392926620982104401176 : [2, 2, 2, 11, 13, 45061, 2338703, 344213870325863]
Time 0.214 sec

完整代码如下：

在线试用！

def FactorECM(N, *, bound = 1 << 8, icurve = 0):
# https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
import math, random
def FactorTrialDiv(x):
# https://en.wikipedia.org/wiki/Trial_division
fs = []
for d in range(2, x + 1):
if d * d > x:
break
while x % d == 0:
fs.append(d)
x //= d
if x != 1:
fs.append(x)
return sorted(fs)
def FermatPRP(n, trials = 32):
# https://en.wikipedia.org/wiki/Fermat_primality_test
for i in range(trials):
if pow(random.randint(2, n - 2), n - 1, n) != 1:
return False
return True
def GenPrimes(end):
# https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
composite = [False] * end
for p in range(2, int(math.sqrt(end) + 3)):
if composite[p]:
continue
for i in range(p * p, end, p):
composite[i] = True
for p in range(2, end):
if not composite[p]:
yield p
def GCD(a, b):
# https://en.wikipedia.org/wiki/Euclidean_algorithm
while b != 0:
a, b = b, a % b
return a
def EGCD(a, b):
# https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
ro, r, so, s = a, b, 1, 0
while r != 0:
ro, (q, r) = r, divmod(ro, r)
so, s = s, so - q * s
return ro, so, (ro - so * a) // b
def PrimePow(p):
bound2 = int(math.sqrt(bound) + 1)
mp = p
while True:
mp *= p
if mp >= bound2:
return mp // p
def ModInv(a, n):
# https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
a %= n
g, s, t = EGCD(a, n)
if g != 1:
raise ValueError(a)
return s % n
def EcAdd(A, B, X0, Y0, X1, Y1):
# https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
if X0 == X1 and Y0 == Y1:
# Double
l = ((3 * X0 ** 2 + A) * ModInv(2 * Y0, N)) % N
x = (l ** 2 - 2 * X0) % N
y = (l * (X0 - x) - Y0) % N
else:
# Add
l = ((Y1 - Y0) * ModInv(X1 - X0, N)) % N
x = (l ** 2 - X0 - X1) % N
y = (l * (X0 - x) - Y0) % N
return x, y
def EcMul(A, B, X, Y, k):
# https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
assert k >= 2, k
k -= 1
BX, BY = X, Y
while k != 0:
if k & 1:
X, Y = EcAdd(A, B, X, Y, BX, BY)
BX, BY = EcAdd(A, B, BX, BY, BX, BY)
k >>= 1
return X, Y
def NextECM(x):
return FactorECM(x, bound = bound + (1 << 8), icurve = icurve + 1)
def Main():
# Main ECM Algorithm    

if N < (1 << 16):
fs = FactorTrialDiv(N)
if len(fs) >= 2:
print('Factors TrialDiv:', fs, flush = True)
return fs

if FermatPRP(N):
return [N]

print(f'Curve {icurve:>4},  bound 2^{math.log2(bound):>6.2f}', flush = True)

while True:
X, Y, A = [random.randrange(N) for i in range(3)]
B = (Y ** 2 - X ** 3 - A * X) % N
if 4 * A ** 3 - 27 * B ** 2 == 0:
continue
break

for p in GenPrimes(bound):
k = PrimePow(p)
try:
X, Y = EcMul(A, B, X, Y, k)
except ValueError as ex:
g = GCD(ex.args[0], N)
assert g > 1, g
if g != N:
print('Factor ECM:', g, flush = True)
return sorted(NextECM(g) + NextECM(N // g))

return NextECM(N)

return Main()
def Test():
import random, time
bits = 72
N = random.randrange(1 << (bits - 1), 1 << bits)
print('N:', N)
tb = time.time()
print(N, ':', FactorECM(N))
print('Time', round(time.time() - tb, 3), 'sec')
if __name__ == '__main__':
Test()