Internet 中存在不同的 Python 分解模块。但是,如果您想自己实现分解(不使用外部库),那么我可以建议您非常快速且非常容易地实现Pollard-Rho Algorithm。我在下面的代码中完全实现了它,如果您不想阅读,您只需直接向下滚动到我的代码(在答案的底部)。
Pollard-Rho 算法很有可能在O(Sqrt(P)) 的时间内找到最小的非平凡因子P(不等于1 或N)。为了比较,您在问题中实现的Trial Division 算法需要O(P) 时间才能找到因子P。这意味着例如,如果一个素数因子 P = 1 000 003 然后试用除法将在 1 000 003 除法运算之后找到它,而 Pollard-Rho 平均会在 1 000 运算之后找到它(Sqrt(1 000 003) = 1 000),这要快得多.
为了使 Pollard-Rho 算法更快,我们应该能够检测素数,将它们排除在因式分解之外,并且不要等待不必要的时间,因为在我的代码中我使用了 Fermat Primality Test,它非常快速且易于只需 7-9 行代码即可实现。
Pollard-Rho算法本身很短,13-15行代码,你可以在我的pollard_rho_factor()函数的最底部看到,剩下的几行代码是辅助函数。
我从头开始实现了所有算法,没有使用额外的库(random 模块除外)。这就是为什么你可以在那里看到我的 gcd() 函数,尽管你可以使用内置 Python 的 math.gcd() 代替(它找到 Greatest Common Divisor)。
您可以在我的代码中看到函数Int(),它仅用于将Python 的整数转换为GMPY2。 GMPY2 int 将使算法更快,您可以使用 Python 的 int(x) 代替。我没有使用任何特定的 GMPY2 函数,只是将所有整数转换为 GMPY2 整数以提高大约 50% 的速度。
例如,我考虑了Pi 的前 190 位数字!!!将它们分解需要 3-15 秒。 Pollard-Rho 算法是随机的,因此每次运行需要不同的时间来分解相同的数字。你可以重新启动程序,看看它会打印出不同的运行时间。
当然,分解时间很大程度上取决于素数除数的大小。一些 50-200 位数字可以在几分之一秒内分解,有些需要几个月。我的示例 190 位 Pi 的素因数非常小,除了最大的一个,这就是它速度快的原因。 Pi 的其他数字可能没有那么快分解。所以数字的位数并不重要,只有素因数的大小很重要。
我有意将pollard_rho_factor() 函数实现为一个独立函数,而不是将其分解为更小的独立函数。虽然它打破了 Python 的风格指南,但(我记得)它建议不要嵌套函数并将所有可能的函数放在全局范围内。样式指南还建议在脚本的第一行中在全局范围内进行所有导入。我特意做了一个函数,以便它可以轻松复制粘贴并完全准备好在您的代码中使用。 Fermat primality test is_fermat_probable_prime() 子函数也是可复制粘贴的,无需额外依赖即可工作。
在极少数情况下,Pollard-Rho 算法可能无法找到非平凡素因子,尤其是对于非常小的因子,例如,您可以将 test() 内的 n 替换为小数字 4 并看到 Pollard-罗失败了。对于如此小的失败因素,您可以轻松使用您在问题中实现的 Trial Division 算法。
Try it online!
def pollard_rho_factor(N, *, trials = 16):
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
import math, random
def Int(x):
import gmpy2
return gmpy2.mpz(x) # int(x)
def is_fermat_probable_prime(n, *, trials = 32):
# https://en.wikipedia.org/wiki/Fermat_primality_test
import random
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 gcd(a, b):
# https://en.wikipedia.org/wiki/Greatest_common_divisor
# https://en.wikipedia.org/wiki/Euclidean_algorithm
while b != 0:
a, b = b, a % b
return a
def found(f, prime):
print(f'Found {("composite", "prime")[prime]} factor, {math.log2(f):>7.03f} bits... {("Pollard-Rho failed to fully factor it!", "")[prime]}')
return f
N = Int(N)
if N <= 1:
return []
if is_fermat_probable_prime(N):
return [found(N, True)]
for j in range(trials):
i, stage, y, x = 0, 2, Int(1), Int(random.randint(1, N - 2))
while True:
r = gcd(N, abs(x - y))
if r != 1:
break
if i == stage:
y = x
stage <<= 1
x = (x * x + 1) % N
i += 1
if r != N:
return sorted(pollard_rho_factor(r) + pollard_rho_factor(N // r))
return [found(N, False)] # Pollard-Rho failed
def test():
import time
# http://www.math.com/tables/constants/pi.htm
# pi = 3.
# 1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679
# 8214808651 3282306647 0938446095 5058223172 5359408128 4811174502 8410270193 8521105559 6446229489 5493038196
# n = first 190 fractional digits of Pi
n = 1415926535_8979323846_2643383279_5028841971_6939937510_5820974944_5923078164_0628620899_8628034825_3421170679_8214808651_3282306647_0938446095_5058223172_5359408128_4811174502_8410270193_8521105559_6446229489
tb = time.time()
print('N:', n)
print('Factors:', pollard_rho_factor(n))
print(f'Time: {time.time() - tb:.03f} sec')
test()
输出:
N: 1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489
Found prime factor, 1.585 bits...
Found prime factor, 6.150 bits...
Found prime factor, 20.020 bits...
Found prime factor, 27.193 bits...
Found prime factor, 28.311 bits...
Found prime factor, 545.087 bits...
Factors: [mpz(3), mpz(71), mpz(1063541), mpz(153422959), mpz(332958319), mpz(122356390229851897378935483485536580757336676443481705501726535578690975860555141829117483263572548187951860901335596150415443615382488933330968669408906073630300473)]
Time: 2.963 sec