使用GMPY2（或GMP）找出所有因子的最高效方法？
我知道已经有类似的问题了，但我想用GMPY2（或者其他类似的GMP工具）来加速这个过程。以下是我现在的代码，效果还不错，但有没有更好的方法呢？
编辑：新代码，检查了2和3这两个因数。
def factors(n):
result = set()
result |= {mpz(1), mpz(n)}
def all_multiples(result, n, factor):
z = mpz(n)
while gmpy2.f_mod(mpz(z), factor) == 0:
z = gmpy2.divexact(z, factor)
result |= {mpz(factor), z}
return result
result = all_multiples(result, n, 2)
result = all_multiples(result, n, 3)
for i in range(1, gmpy2.isqrt(n) + 1, 6):
i1 = mpz(i) + 1
i2 = mpz(i) + 5
div1, mod1 = gmpy2.f_divmod(n, i1)
div2, mod2 = gmpy2.f_divmod(n, i2)
if mod1 == 0:
result |= {i1, div1}
if mod2 == 0:
result |= {i2, div2}
return result
如果可能的话，我也想知道有没有只在n^(1/3) 和 2^(2/3)*n(1/3)范围内检查因数的实现。
举个例子，Mathematica的factor()函数比我写的Python代码快得多。我想对20到50位的十进制数字进行因数分解。我知道ggnfs可以在不到5秒的时间内完成这个任务。
我想了解一下，Python中是否也有实现快速因数分解的模块。
2 个回答
网上有很多不同的Python因式分解模块。如果你想自己实现因式分解（不使用外部库），我可以推荐一个相对快速且容易实现的算法——Pollard-Rho算法。如果你不想阅读，可以直接滚动到我下面的代码部分（在回答的底部）。
Pollard-Rho算法很有可能在时间复杂度为O(Sqrt(P))的情况下找到最小的非平凡因子P（不等于1或N）。相比之下，你在问题中实现的试除法算法需要O(P)的时间来找到因子P。举个例子，如果一个质因子P = 1 000 003，那么试除法需要进行1 000 003次除法操作才能找到它，而Pollard-Rho算法平均只需1 000次操作（Sqrt(1 000 003) = 1 000），速度快得多。
为了让Pollard-Rho算法更快，我们需要能够检测质数，以排除它们的因式分解，避免不必要的等待。为此，我在代码中使用了费马质数测试，这个测试非常快，且只需7-9行代码就能实现。
Pollard-Rho算法本身代码很短，大约13-15行，你可以在我pollard_rho_factor()函数的最底部看到，剩下的代码是一些辅助函数。
我从头开始实现了所有算法，没有使用额外的库（除了random模块）。所以你可以看到我的gcd()函数，尽管你也可以使用Python内置的math.gcd()函数（用于找到最大公约数）。
你可以在我的代码中看到Int()函数，它只是用来将Python的整数转换为GMPY2。使用GMPY2的整数会让算法更快，你也可以直接使用Python的int(x)。我没有使用任何特定的GMPY2函数，只是将所有整数转换为GMPY2的整数，这样可以提高大约50%的速度。
举个例子，我对圆周率的前190位数字进行了因式分解！这需要3到15秒的时间。由于Pollard-Rho算法是随机的，因此每次运行同一个数字的因式分解所需的时间可能不同。你可以重新启动程序，看看它会打印出不同的运行时间。
当然，因式分解的时间很大程度上取决于质因子的大小。一些50-200位的数字可以在瞬间完成因式分解，而有些则可能需要几个月。我的例子中，圆周率的190位数字有相对较小的质因子，除了最大的一个，因此分解速度很快。圆周率的其他数字可能就没有那么快了。所以，数字的位数并不是最重要的，质因子的大小才是关键。
我故意将pollard_rho_factor()函数实现为一个独立的函数，没有将其拆分成更小的函数。虽然这样违反了Python的风格指南（我记得它建议不要有嵌套函数，并将所有可能的函数放在全局作用域中），而且风格指南还建议在脚本的前几行进行所有的导入。我这样做是为了让它更容易复制粘贴，并且可以直接在你的代码中使用。费马质数测试的is_fermat_probable_prime()子函数也可以直接复制粘贴，并且不需要额外的依赖。
在极少数情况下，Pollard-Rho算法可能无法找到非平凡的质因子，尤其是对于非常小的因子。例如，你可以将test()中的n替换为小数字4，看看Pollard-Rho会失败。对于这样的小因子，你可以轻松使用你在问题中实现的试除法。
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
我对你的代码做了一些快速修改，去掉了一些重复的名字查找。算法还是那个算法，但在我的电脑上运行速度快了大约一倍。
import gmpy2
from gmpy2 import mpz
def factors(n):
result = set()
n = mpz(n)
for i in range(1, gmpy2.isqrt(n) + 1):
div, mod = divmod(n, i)
if not mod:
result |= {mpz(i), div}
return result
print(factors(12345678901234567))
其他建议需要更多关于数字大小的信息等等。例如，如果你需要所有可能的因子，从所有的质因子构造这些因子可能会更快。这个方法可以让你在进行时减少范围声明的限制，并且在去掉所有的2的因子后，你还可以每次加2。
更新 1
我对你的代码做了一些额外的修改。我觉得你的 all_multiplies() 函数不太对。你的 range() 声明不是最优的，因为2又被检查了一次，而我第一次修复反而让它变得更糟。
新代码会等到知道余数是0时再计算共因子。我还尽量多用内置函数。例如，mpz % integer 比 gmpy2.f_mod(mpz, integer) 或 gmpy2.f_mod(integer, mpz) 要快，其中 integer 是普通的Python整数。
import gmpy2
from gmpy2 import mpz, isqrt
def factors(n):
n = mpz(n)
result = set()
result |= {mpz(1), n}
def all_multiples(result, n, factor):
z = n
f = mpz(factor)
while z % f == 0:
result |= {f, z // f}
f += factor
return result
result = all_multiples(result, n, 2)
result = all_multiples(result, n, 3)
for i in range(1, isqrt(n) + 1, 6):
i1 = i + 1
i2 = i + 5
if not n % i1:
result |= {mpz(i1), n // i1}
if not n % i2:
result |= {mpz(i2), n // i2}
return result
print(factors(12345678901234567))
我建议你把程序改成只找出小于 n 平方根的所有质因子，然后再构造所有的共因子。每次找到一个因子后，就减少 n 的值，检查 n 是否是质数，只有在 n 不是质数的情况下才继续寻找更多因子。
更新 2
pyecm 模块应该能够处理你想要分解的大小数字。下面的例子大约在一秒内完成。
>>> import pyecm
>>> list(pyecm.factors(12345678901234567890123456789012345678901, False, True, 10, 1))
[mpz(29), mpz(43), mpz(43), mpz(55202177), mpz(2928109491677), mpz(1424415039563189)]