感谢您提出非常有趣的问题！
我对您的代码进行了一些优化，将其速度提高了大约 95 倍 倍，在它们对 166 位数字（50 个十进制数字，如您所愿）进行因式分解后需要 10 分钟（请参阅 @987654321 @) 到我的笔记本电脑上的 4 小时（时间取决于运气）和 132 位（40 个十进制数字）只需要 1 到 15 分钟。请参阅本文末尾的最终代码。
进行了以下改进：
-
使用 multiprocessing 模块占用所有 CPU 内核，而不是您的代码的单核版本。这项更改使您的代码在我的 4 核笔记本电脑（具有 8 个硬件线程）上的速度提高了 8 倍。
-
使用GMPY2 库将纯Python int 类型替换为gmpy2.mpz() 类型，由于高度优化的GMP 库，它的所有数学运算速度比Python 快2-3 倍。可以通过注释行 import gmpy2 来禁用 gmpy2 模块的使用（如果您愿意）。
-
将纯 Python inv() 函数的主体替换为基于 GMPY2 的快速 gmpy2.gcdext() 实现 Extended Euclidean Algorithm。如果在我的代码中禁用了 gmpy2 使用（这是可能的），则内置 Python 的 pow(a, -1, n) 用于反转整数。 gmpy2.gcdext() 比 Python 的 pow(a, -1, n) 快，并且比您的扩展欧几里得算法的 inv() 实现快得多。在 Python 的 C 源代码中，通过 pow(a, -1, n) 反转数字的方法与您的扩展欧几里得算法类似，您可以使用 see here 的源代码，这个版本比您的版本快，唯一的原因是它完全基于 C，原因GMPY2 的变体gmpy2.gcdext() 速度更快，因为它使用了更深奥的数学算法。
-
将曲线点乘以range(2, bound) 替换为仅乘以[2 ; bound] 范围内的所有素数。这是一个显着的改进，因为素数比所有数字的范围要稀疏得多，因此乘法次数要少得多。为了使这种改变成为可能，我必须引入bound_pow 常量，使得prime ** m <= bound_pow，换句话说，我乘以它的不是素数本身，而是乘以它在一定范围内的幂。当曲线的阶数大于某些素因子的 1 次方时，需要用到素数的幂。
-
删除了不必要的inv(P.x, n) 和inv(P.y, n)，它们占用了大量时间但没有做任何有用的工作，因为这两个反转产生因子的机会并不比反转任何随机数大，所以根本没有用。
-
删除了 ECpoint 构造函数主体中对椭圆曲线不变量 (y ** 2 - x ** 3 - A * x - B) % N == 0 的不必要检查，因为此检查在每次 __add__() 调用上完成，并且仅当点添加执行不正确时才会失败，因此它只是一个减慢速度的调试检查放下东西。
-
做了很多可视化，打印到不同值、时间和进度表的控制台。
-
实现了用于计算许多有用帮助值的函数，例如 bound 的最佳值，而不是您的近似公式（我没有使用，只是打印了）。还计算了所需曲线的最佳数量。还用于计算给定边界和曲线数量的最大因子大小（错过10% 的概率）。计算给定边界和曲线数量所需的工作量 (Work2()) 的函数。
上面的步骤 1)-6) 都在提高速度（步骤 5)-6) 不是很重要，只有 1)-4) 很重要），您可以在代码中手动重现它们，它们非常简单.
112 位复合的示例输出（在 148 秒内分解）：
Factoring 2524582745259710504267101693459139
Calculating. Wait...
n 2^110.96, bound 2^13.18 (optimal 2^14.24, 1.051x faster),
need at most 974 curves (factor up to 2^55.48)
ooo
curves 1/974 ( 0.103%), factors < 2^14.508, time 00:00:20
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
curves 86/974 ( 8.830%), factors < 2^41.646, time 00:01:20
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
curves 179/974 ( 18.378%), factors < 2^45.995, time 00:02:21
oooooo
Factored on curve 190 (19.61%) (rand_seed = 248080385138786)!
ooo
Factors: P 38240987134592087, P 66017719060814197
2 prime, 0 composite
All factors are prime!!!
Time 148.2 sec
166 位数字（50 位）的示例，使用 random.seed(1)，需要 11.5 分钟来计算：
Factoring 45733304697523851846830687775886905451041736136239
Calculating. Wait...
n 2^164.97, bound 2^16.79 (optimal 2^17.43, 1.023x faster), need at most 5722 curves (factor up to 2^82.48)
........................................o
curves 1/5722 ( 0.017%), factors < 2^18.480, time 00:00:25
ooooooo........................................oooooooo........................................oooooooo........................................o
curves 25/5722 ( 0.437%), factors < 2^43.197, time 00:01:25
ooooooo........................................oooooooo........................................oooooooo......................................o
curves 49/5722 ( 0.856%), factors < 2^48.646, time 00:02:26
..ooooooo....................................o....ooooooo....................................o....ooooooo....................................o
curves 73/5722 ( 1.276%), factors < 2^51.779, time 00:03:27
...o.oooooo..................................o.....o.oooooo..................................o....o.o.ooooo..................................o
curves 97/5722 ( 1.695%), factors < 2^53.972, time 00:04:28
....o.oo.oooo..................................o....o.oooo.oo..................................o....o.ooo.ooo..................................o
curves 121/5722 ( 2.115%), factors < 2^55.655, time 00:05:29
.....oooo.ooo..................................o....o.ooo.oo.o.................................o....o.ooo.oo.o.................................o
curves 145/5722 ( 2.534%), factors < 2^57.017, time 00:06:29
....o.oooo.o.o.................................o....o.ooo.oo.o.................................o...o..oooo.o.o.................................o
curves 169/5722 ( 2.954%), factors < 2^58.161, time 00:07:30
o.....oooo.o.o.................................oo.....oooo.o.o.................................oo.....oo.ooo.o.................................o
curves 193/5722 ( 3.373%), factors < 2^59.145, time 00:08:31
o.....oo.ooo.o.................................oo.....oo.ooo.o.................................oo.....oooo..oo................................o
curves 217/5722 ( 3.792%), factors < 2^60.009, time 00:09:31
.o.....oooo..oo.................................oo.....oooo..o.o................................oo.....oooo..o.o...............................o
curves 241/5722 ( 4.212%), factors < 2^60.778, time 00:10:32
.o.....oooo..o.o...............................o.o.....oooo..o.o...............................o.o.....oooo..o.o...............................o
curves 265/5722 ( 4.631%), factors < 2^61.470, time 00:11:33
.o.....oooo..o.o...............................o.o....o.ooo..o.o............
Factored on curve 282 (4.95%) (rand_seed = 83739500771956)!
Factors: P 4856082635476451488745861, P 9417736091106851081336099
2 prime, 0 composite
All factors are prime!!!
Time 724.9 sec
为了在我的test() 函数中进行演示，我将一个复合数组合为两个位长相等的素数的乘积（与 RSA 中相同）。合数的位长由bits常量给出。
由于在脚本开头使用了random.seed(0)，我的代码中的所有随机性都是确定性的。此种子已设置为固定，以便您可以在多次运行脚本中重现完全相同的确定性结果。如果您希望有其他随机性并因此产生其他合成和曲线，您可以更改此种子值。正如您在上面的示例控制台输出中看到的那样，最后我显示了对一个数字进行分解的成功曲线的rand_seed 值，如果您想再次重现曲线的精确参数以快速重做分解，则此种子值很有用，或者把这条成功的曲线送给你的朋友。
计算所需的曲线数，如果一个因子的位长等于复合位长的一半，那么您恰好有10% 的机会错过一个因子。换句话说，在最坏的情况下，当一个因子可能最大时，如果你运行所有这些曲线数量，那么你恰好有10% 错过这个因子的机会（90% 成功分解的机会）。如果您再跑两次曲线，那么您将获得1% 错过机会（99% 成功机会）。如果您想获得50% 分解的机会，那么您只需检查1/3-th 的计算曲线数量。很多时候，我在检查1/10-th 所需的曲线数量后就成功分解了，在这种情况下，达到一个因子的机会是10%（90% 错过机会）。
与纯 Python 整数相比，使用 GMPY2 库提供了大约 1.8x 倍的整体加速，非常显着，几乎快了两倍。如果您想要纯 Python 解决方案，那么只需在脚本的第一行注释掉 import gmpy2。在 Linux 上，您可以通过 sudo apt install python3-gmpy2 安装 gmpy2。在 Windows 上下载相应 Python 版本的轮子文件 from here 并通过例如安装它python -m pip install gmpy2‑2.0.8‑cp39‑cp39‑win_amd64.whl.
如果您计算Non-Adjacent Form, wNAF，则可以进一步提高20% 的速度，这种方法大大减少了所需的加法量（从50% 用于1 位的随机分布到10%），加倍的数量（添加相同的指向同一点）保持不变。我在我的ECM 库（其他脚本）中使用了这个方法，它给了我20% 的整体加速。此外，要使用此方法，您必须连续预先计算（乘）围绕200 素数的乘积，这与当前只能将素数相乘的方法不同。
您可以看到我预先计算并打印到控制台optimal 绑定，使用OptimalBoundLog() 函数计算，使用此绑定应该给3-5% 更多的加速，尽管我没有在下面的代码中使用此绑定。
注意。我还在 Python 中实现了我自己的椭圆曲线分解脚本，该脚本相当长，并且它在纯 C++（使用 Cython 包装器）中实现了点乘法。我做了纯 Python（基于 GMPY2）乘法和 C++ 版本的计时，它们不同 3x 次。所以 C++ 提供了额外的3x 倍提升。
如果要近似计算总加速比，那么我们将有：8x 由于多处理（我家里有 8 个内核），1.8x 由于使用 GMPY2，2x 由于使用 gmpy2.gcdext()而不是 Pythonic 版本的扩展欧几里得算法，在 3x 附近，因为使用素数而不是连续数字，1.1x 因为删除了所有不必要的代码，如 inv(P.x, n)，3x，因为使用 C++ 进行点乘法（我只在我的库中使用过，这里没有发布），如果实现Non-Adjacent Form (wNAF)，也可以获得1.2x boost。总共 340 倍 倍的加速！
您还应该考虑到实际运行时间很大程度上取决于您的运气。有时你应该在找到一个因素之前检查所有计算出的所需曲线数量的3%（如果你很幸运的话），有时你必须检查它们中的200%（如果你很不走运的话）。检查其中的100% 将给出90% 找到一个因素的概率（或10% 错过一个因素的概率）。当然时间也取决于一个最小因素的比特长度。所需曲线的数量是针对最坏情况计算的，即因子的位长等于输入合数n 的位长的一半。
您可以在成功分解的最后看到我显示类似于Factored on curve 190 (19.61%) (rand_seed = 248080385138786)! 的内容，这里rand_seed 非常重要，它唯一标识了成功分解该数字的最后一条曲线。如果您长时间运行分解，那么最后您应该记住这个种子值，以便您可以立即重现分解的结果。要重现结果，只需在 ECM() 函数中输入相同的复合数，并将代码中的 'rand_seed': random.randrange(1 << 48), 替换为 'rand_seed': 248080385138786,（此种子值取自上面的示例输出）。然后重新运行脚本，您应该会在第 0 条曲线的最开始看到成功的因式分解结果。非常重要 - 仅暂时在脚本中对rand_seed 进行此修改，完成后将其恢复为'rand_seed': random.randrange(1 << 48),，否则您将破坏其他输入@987654409 的所有随机性来源@ 数字，这将导致无限次不成功的因式分解运行。您也可以将此rand_seed 作为ECM() 函数的可选参数，这样您就无需就地修改代码。
要使用以下代码尝试修改两个值 - 输入复合数的位长度，它是通过 test() 函数内脚本末尾的 bits 变量设置的（现在是 112 位）。第二 - 尝试在脚本的最开始修改行random.seed(0)，将种子值更改为其他值，如1、2、3 等。如果你不更改这个种子，那么你会得到完全每次运行脚本的结果相同。此种子控制脚本内所有随机值的行为。如果您出于某种原因想要重现完全相同的脚本运行结果（例如向您的朋友显示结果），则需要相同的种子。
改进的代码本身：
Try it online!
gmpy2 = None
import gmpy2
import random
random.seed(0)
class InvError(Exception):
def __init__(self, v):
self.value = v
def Int(x):
return int(x) if gmpy2 is None else gmpy2.mpz(x)
def inv(a, n):
a %= n
if gmpy2 is None:
try:
return pow(a, -1, n)
except ValueError:
import math
raise InvError(math.gcd(a, n))
else:
g, s, t = gmpy2.gcdext(a, n)
if g != 1:
raise InvError(g)
return s % n
# Slow version below
'''
r1, s1, t1 = 1, 0, a
r2, s2, t2 = 0, 1, n
while t2:
q = t1//t2
r1, r2 = r2, r1-q*r2
s1, s2 = s2, s1-q*s2
t1, t2 = t2, t1-q*t2
if t1!=1: raise InvError(t1)
else: return r1
'''
class ECpoint(object):
def __init__(self, A, B, N, x, y, *, prepare = True):
if prepare:
N = Int(N)
A, B, x, y = [Int(e) % N for e in [A, B, x, y]]
if (y ** 2 - x ** 3 - A * x - B) % N != 0:
raise ValueError
self.A, self.B, self.N, self.x, self.y = A, B, N, x, y
def __add__(self, other):
A, B, N = self.A, self.B, self.N
Px, Py, Qx, Qy = self.x, self.y, other.x, other.y
if Px == Qx and Py == Qy:
s = ((Px * Px * 3 + A) * inv(Py * 2, N)) % N
else:
s = ((Py - Qy) * inv(Px - Qx, N)) % N
x = (s * s - Px - Qx) % N
y = (s * (Px - x) - Py) % N
return ECpoint(A, B, N, x, y, prepare = False)
def __rmul__(self, other):
other = Int(other - 1)
r = self
while True:
if other & 1:
r = r + self
if other == 1:
return r
other >>= 1
self = self + self
def BinarySearch(f, a, b):
while a < b:
m = (a + b) // 2
if f(m):
b = m
else:
a = m + 1
assert a == b and f(a), (a, b, f(a))
return a
def FactorBitSize(bound, ncurves, miss_prob = 0.1):
import math
bound_log2 = math.log2(bound)
def Prob(factor_log2):
x = factor_log2 / bound_log2
return x ** -x
def F(factor_log2):
return (1. - Prob(factor_log2)) ** ncurves >= miss_prob
return BinarySearch(lambda x: F(x / 1000.), math.log2(bound) * 1000., 512 * 1000.) / 1000.
def NeededCurves(bound, target_fac_log2, miss_prob = 0.1):
def F(ncurves):
return FactorBitSize(bound, ncurves, miss_prob) >= target_fac_log2
return round(BinarySearch(lambda x: F(x), 1, 10 ** 15))
def Work(bound, bound_pow, *, logs = [0.], cache = {}):
if bound not in cache:
import math
for ilog in range(100):
if get_prime(1 << ilog) >= bound:
break
cnt_primes = BinarySearch(lambda i: get_prime(i) >= bound, 0, 1 << ilog)
bound_pow = max(bound, bound_pow)
bound_log = math.log2(bound)
bound_pow_log = math.log2(bound_pow)
while cnt_primes >= len(logs):
plog = math.log2(get_prime(len(logs)))
sum_plog = plog
while True:
sum_plog += plog
if sum_plog >= max(bound_log, bound_pow_log):
sum_plog -= plog
break
logs.append(logs[-1] + sum_plog)
cache[bound] = logs[cnt_primes]
return cache[bound]
def Work2(bound, bound_pow, factor_log):
return Work(bound, bound_pow) * NeededCurves(bound, factor_log)
def OptimalBoundLog(factor_log, bound_pow):
import math
mwork = None
bound_log_start = 10
for bound_log in range(bound_log_start, math.floor(factor_log) + 1):
bound = 2 ** bound_log
work = Work(bound, bound_pow) * NeededCurves(bound, factor_log)
if mwork is not None and work > mwork[0]:
break
if mwork is None or work < mwork[0]:
mwork = (work, bound_log)
else:
bound_log = bound_log_start + 1
mult = 200
mwork = None
for bound_log2 in range((bound_log - 3) * mult, bound_log * mult + 1):
bound = round(2 ** (bound_log2 / mult))
work = Work(bound, bound_pow) * NeededCurves(bound, factor_log)
if mwork is None or work < mwork[0]:
mwork = (work, bound_log2 / mult)
return mwork[1]
def get_prime(i, *, primes = [2, 3]):
while i >= len(primes):
for n in range(primes[-1] + 2, 1 << 62, 2):
isp = True
for p in primes:
if p * p > n:
break
if n % p == 0:
isp = False
break
if isp:
primes.append(n)
break
return primes[i]
def prime_power(idx, bound, bound_pow, *, cache = {}):
key = (bound, bound_pow)
if key not in cache:
bound_pow = max(bound, bound_pow)
r = []
for i in range(1 << 62):
p = get_prime(i)
if p >= bound:
break
m = p
while True:
m2 = m * p
if m2 >= bound_pow:
break
m = m2
r.append(m)
cache[key] = r
return cache[key][idx] if idx < len(cache[key]) else None
def ExcInfo(ex):
return f'{type(ex).__name__}: {ex}'
def ProcessCurve(*, n, bound, bound_pow, shared, rand_seed, curve_idx, num_curves):
try:
random.seed(rand_seed)
x0, y0, a = [random.randrange(1, n) for i in range(3)]
# x0 = 2; y0 = 3
b = (y0 ** 2 - x0 ** 3 - a * x0) % n
P = ECpoint(a,b,n, x0,y0)
for i in range(1 << 62):
if shared.get('finish', False):
return {'ex': 'Interrupted: Finishing...'}
k = prime_power(i, bound, bound_pow)
if i > 0 and i % 2000 == 0 or k is None:
Print(('.', 'o')[k is None], end = '', flush = True)
if k is None:
break
P = k * P
except InvError as e:
d = e.value
if d != n:
Print(f'\nFactored on curve {curve_idx} ({(curve_idx + 1) / num_curves * 100.:.02f}%) (rand_seed = {rand_seed})!')
return {'factors': sorted([d, n // d])}
except BaseException as ex:
shared['finish'] = True
return {'ex': ExcInfo(ex)}
return {}
def Print(*pargs, **nargs):
print(*pargs, **{'flush': True, **nargs})
def ECM(n, *, processes = None):
import math, multiprocessing as mp, time
from datetime import datetime as dt
Print(f'Factoring {n}')
if fermat_prp(n):
return [n]
Print('Calculating. Wait...')
start_time = dt.now()
bound = max(int(math.e**(1/2*math.sqrt(math.log(n)*math.log(math.log(n))))),100)
bound_pow = max(bound, 1 << 18)
factor_log = math.log2(n) / 2
max_curves = NeededCurves(bound, factor_log)
opt_bound_log = OptimalBoundLog(factor_log, bound_pow)
processes = processes or mp.cpu_count()
Print(f'n 2^{math.log2(n):.02f}, bound 2^{math.log2(bound):.02f} (optimal 2^{opt_bound_log:.02f}, ' +
f'{Work2(bound, bound_pow, factor_log) / Work2(2 ** opt_bound_log, bound_pow, factor_log):.03f}x faster), ' +
f'need at most {max_curves} curves (factor up to 2^{factor_log:.02f})')
with mp.Manager() as manager, mp.Pool(processes) as pool:
try:
ncurves, report_time = 0, 0
shared = manager.dict()
res = []
for icurve in range(1 << 62):
res.append(pool.apply_async(ProcessCurve, (),
{
'n': n, 'bound': bound, 'bound_pow': bound_pow,
'shared': shared, 'rand_seed': random.randrange(1 << 48),
'curve_idx': icurve, 'num_curves': max_curves,
}))
if len(res) < processes * 9:
continue
while len(res) >= processes * 6:
res2 = []
for e in res:
if not e.ready():
res2.append(e)
continue
e = e.get()
assert 'ex' not in e, e['ex']
if 'factors' in e:
return e['factors']
ncurves += 1
if time.time() - report_time >= 60:
Print(f'\ncurves {ncurves:>5}/{max_curves} ({ncurves / max_curves * 100.:>7.03f}%), factors ' +
f'< 2^{FactorBitSize(bound, ncurves):.03f}, ' +
f'time {(dt(2000, 1, 1) + (dt.now() - start_time)).strftime("%H:%M:%S"):>9}')
report_time = time.time()
res = res2
time.sleep(0.01)
except BaseException as ex:
Print(f'\nException: {ExcInfo(ex)}. Finishing, wait!')
finally:
shared['finish'] = True
pool.close()
pool.join()
Print()
return [n]
def fermat_prp(n, trials = 32):
# https://en.wikipedia.org/wiki/Fermat_primality_test
if n <= 16:
return n in (2, 3, 5, 7, 11, 13)
for i in range(trials):
if pow(random.randint(2, n - 2), n - 1, n) != 1:
return False
return True
def gen_random_prime(bits):
while True:
n = random.randrange(1 << (bits - 1), 1 << bits)
if fermat_prp(n):
return n
def Prod(it):
import functools
return functools.reduce(lambda x, y: x * y, it, 1)
def test():
import time
# First 190 digits of Pi
pi = 1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489
bits = 112
nprimes = 2
n = Prod([gen_random_prime(bits // nprimes) for i in range(nprimes)])
tb = time.time()
fs = ECM(n)
Print('Factors:', ', '.join([('C', 'P')[fermat_prp(e)] + f' {e}' for e in fs]))
assert n == Prod(fs), (n, fs)
Print(sum(fermat_prp(e) for e in fs), 'prime,', sum(not fermat_prp(e) for e in fs), 'composite')
Print('All factors are prime!!!' if all(fermat_prp(e) for e in fs) else 'Composite factors remaining...')
Print(f'Time {time.time() - tb:.01f} sec')
if __name__ == '__main__':
test()