如何计算 Python 中的非常长的数字
How to calculate very long numbers in Python
我正在尝试使用 python 计算 Pollard rho 数字以获得非常长的整数，例如小于一
65779646778470582047547160396995720887221575959770627441205850493179860146690755880473849736210807898494458426111244201404810495587574110361900128405354081638434164434968839614760264675889940272767106444249
我已经尝试计算我的 intel core i9 10980HK CPU，这导致几分钟的高负载工作没有任何成功。我正在尝试使用 numba 和 @njit 装饰器来连接 RTX 2070 super（在笔记本电脑上），但它给出了以下错误。
- argument 0: Int value is too large:
这里是代码：
import numpy as np
import datetime
def pgcd(a,b):
if b==0:
return a
else:
r=a%b
return pgcd(b,r)
def pollardrho(n):
f = lambda z: z*z+1
x, y, d = 1, 1, 1
c = 0
while d==1:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = pgcd(y-x, n)
return d, c
def test_time(n):
t = datetime.datetime.now()
d, c = pollardrho(int(n))
tps = datetime.datetime.now() - t
print(tps, c, d)
file = open("some_powersmooths_large.txt", "r")
for line in file.readlines():
if not line.startswith("#"):
print(line)
print("\n")
test_time(line)
print("\n")
我该如何处理这种大数计算。
鉴于pollardrho中的操作效率很低，操作需要一段时间我并不奇怪。
但是，我不知道具体的功能，所以我不知道是否可以提高效率。
在Python中，整数具有任意长度。
这意味着它们可以是任意长度，并且 Python 本身将处理使用 64 位整数正确存储它（通过将它们分布在多个整数上）。
（你可以自己测试这个，例如创建一个不能存储在 64 位无符号整数中的整数，比如 a = 2**64，然后检查 a.bit_length() 方法的输出，它应该是 65)
因此，从理论上讲，您应该能够计算任何整数。
但是，由于您使用的是 Numba，因此由于 Numba 的工作方式，您只能使用实际上可以存储在 64 位无符号整数中的整数。
您得到的错误只是数字变得太大而无法存储在 64 位无符号整数中。
底线：没有 Numba，您可以很好地计算出该数字。使用 Numba，你不能。
当然，如果你只想粗略地知道数字是多少，而不是精确地知道，你可以直接使用浮点数。
第 1 部分（共 2 部分，请参阅下面的第 2 部分）。
Numba 最多只能处理 64 位整数，它没有大整数运算，只有 Python 有。 Numba 的开发人员承诺，未来版本将支持大整数。您需要大整数运算，因为您的输入和计算中有非常大的整数。
您的一个优化建议是使用 GMPY2 Python library. It is highly-optimized library of long arithmetics, considerably faster than regular Python implementation of long arithmetics. For very large integers for example it implements multiplication using Fast Fourier Transform，这是目前最快的乘法算法。
但是 GMPY2 安装起来有点困难。 Windows 的最新预编译版本可用 by this link。下载适用于您的 Python 版本的 .whl 文件并通过 pip 安装它，例如对于我的 Windows 64 位 Python 3.7，我下载并安装了 pip install gmpy2-2.0.8-cp37-cp37m-win_amd64.whl。对于 Linux，通过 sudo apt install -y python3-gmpy2.
安装最简单
使用 GMPY2 后，您的代码将变得尽可能快，因为这个库代码几乎是世界上最快的。即使是 Numba（如果它有长算术）也不会改进更多。只有更快的公式和更好的算法才能帮助进一步改进或更小的输入整数。
但是即使使用 GMPY2，您的示例大整数对于您的算法来说也是一种变大的方法。您必须选择更小的整数或更快的算法。我已经 运行 你的算法和数字 5 分钟或更长时间，但没有得到结果。但是，如果使用常规 Python 之前的结果会在 1 小时内完成，那么在使用 GMPY2 之后，它可能会在 10 分钟或更快的时间内完成。
也不是很确定，但可能在你的算法中 f(f(y)) % n 应该等同于 f(f(y) % n) % n 应该计算得更快，因为它会做两次较短的乘法。但这需要额外检查。
此外，您的大整数似乎是质数，正如 Primo elliptic curve based primality proving program, it proved primality of this integer in 3 seconds on my PC. Primo only proves primality (with 100% guarantee) but doesn't factor the number (splitting into divisors). Factoring numbers can be done by programs from this list 所证明的那样，这些程序实现了已知最快的因式分解算法，如果某些链接失效，则 Google 那些程序的名称。
只需将所有整数 n 包装成 gmpy2.mpz(n)。例如，我稍微改进了您的代码，将其包装到 gmpy2.mpz() 中，还制作了一个循环，以便打印所有除数。另外作为一个例子，我没有使用你的大质数，而是使用了一个更小的 - Pi 的前 25 位数字，它是复合的，它的所有除数都在我的 PC 上在 7 秒内打印：
import datetime, numpy as np, gmpy2
def num(n):
return gmpy2.mpz(n)
zero, one = num(0), num(1)
def pgcd(a, b):
if b == zero:
return a
else:
r = a % b
return pgcd(b, r)
def pollardrho(n):
f = lambda z: z * z + one
x, y, d = one, one, one
c = 0
while d == one:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = pgcd(y - x, n)
return d, c
def test_time(n):
n = num(int(n))
divs = []
while n > 1:
t = datetime.datetime.now()
d, c = pollardrho(num(int(n)))
tps = datetime.datetime.now() - t
print(tps, c, d, flush = True)
divs.append(d)
assert n % d == 0, (n, d)
n = n // d
print('All divisors:\n', ' '.join(map(str, divs)), sep = '')
test_time(1415926535897932384626433)
#test_time(65779646778470582047547160396995720887221575959770627441205850493179860146690755880473849736210807898494458426111244201404810495587574110361900128405354081638434164434968839614760264675889940272767106444249)
输出：
0:00:00 2 7
0:00:00 10 223
0:00:00.000994 65 10739
0:00:00.001999 132 180473
0:00:07.278999 579682 468017117899
All divisors:
7 223 10739 180473 468017117899
第 2 部分
阅读维基百科文章 (here and here)，我决定实施更快版本的 Pollard-Rho 算法。
我在下面实现的版本看起来更复杂，但除法和乘法少了两倍，平均而言，循环的总迭代次数也少了。
此改进导致我的测试用例 运行ning 时间 3 分钟，与原始 OP 算法 运行ning 时间 7 分钟相比，我的笔记本电脑。
我的算法也是随机的，这意味着它采用 [1, N - 2] 范围内的随机见证，而不是见证 1 或 2。在极少数情况下，它可能会像维基百科中所说的那样失败，然后我用不同的证人重新 运行 算法。它还使用 Fermat primality test 检查输入数字是否为质数，然后不再搜索任何除数。
为了测试，我使用了由代码 p = 1; for i in range(256): p *= random.randrange(2, 1 << 32) 生成的输入数字 p，基本上它由 256 个因素组成，每个因素最多 32-bits。
我还改进了这两种算法以输出更多统计信息。统计参数之一是 pow，它显示了每个步骤的复杂性，0.25 的 pow 表示如果除数是 d，则当前因式分解步骤花费 c = d^0.25 次迭代来找到该除数 d.正如维基百科中所述，Pollard-Rho 算法平均应具有 pow = 0.25，这意味着找到任何除数 d 的复杂度（迭代次数）约为 d^0.25。
在接下来的代码中还有其他改进，例如提供大量统计信息。并找到循环中的所有因素。
我的测试用例算法版本的平均 pow 为 0.24，原来的先前版本有 0.3。较小的 pow 意味着平均进行较少的循环迭代。
还测试了使用和不使用 GMPY2 的我的版本。显然 GMPY2 对常规 Python 大整数算法没有太大改进，主要是因为 GMPY2 针对真正的大数字（数万位）（使用快速傅立叶变换乘法等）进行了更多优化，而这里的数字并不大在我的测试中。但是 GMPY2 仍然提供了大约 1.35x 倍的加速，提供了 运行 宁时间 3 分钟，而没有 GMPY2 的相同算法几乎 4 分钟。要使用或不使用 gmpy2 进行测试，您只需将 def num(n) 函数内部更改为 return gmpy2.mpz(n) 或 return n.
import datetime, numpy as np, gmpy2, random, math
random.seed(0)
def num(n):
return gmpy2.mpz(n)
zero, one = num(0), num(1)
def gcd(a, b):
while b != zero:
a, b = b, a % b
return a
def pollard_rho_v0(n):
f = lambda z: z * z + one
n, x, y, d, c, t = num(n), one, one, one, 0, datetime.datetime.now()
while d == one:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = gcd(y - x, n)
return d, {'i': c, 'n_bits': n.bit_length(), 'd_bits': round(math.log(d) / math.log(2), 2),
'pow': round(math.log(max(c, 1)) / math.log(d), 4), 'time': str(datetime.datetime.now() - t)}
def is_fermat_prp(n, trials = 32):
n = num(n)
for i in range(trials):
a = num((3, 5, 7)[i] if i < 3 else random.randint(2, n - 2))
if pow(a, n - 1, n) != 1:
return False
return True
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
# https://ru.wikipedia.org/wiki/%D0%A0%D0%BE-%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%9F%D0%BE%D0%BB%D0%BB%D0%B0%D1%80%D0%B4%D0%B0
def pollard_rho_v1(N):
AbsD = lambda a, b: a - b if a >= b else b - a
N, fermat_prp, t = num(N), None, datetime.datetime.now()
SecsPassed = lambda: (datetime.datetime.now() - t).total_seconds()
for j in range(8):
i, stage, y, x = 0, 2, num(1), num(random.randint(1, N - 2))
while True:
if (i & 0x3FF) == 0 and fermat_prp is None and (SecsPassed() >= 15 or j > 0):
fermat_prp = is_fermat_prp(N)
if fermat_prp:
r = N
break
r = gcd(N, AbsD(x, y))
if r != one:
break
if i == stage:
y = x
stage <<= one
x = (x * x + one) % N
i += 1
if r != N or fermat_prp:
return r, {'i': i, 'j': j, 'n_bits': N.bit_length(), 'd_bits': round(math.log(r) / math.log(2), 2),
'pow': round(math.log(max(i, 1)) / math.log(r), 4), 'fermat_prp': fermat_prp, 'time': str(datetime.datetime.now() - t)}
assert False, f'Pollard-Rho failed after {j + 1} trials! N = {N}'
def factor(n, *, ver = 1):
assert n > 0, n
n, divs, pows, tt = int(n), [], 0., datetime.datetime.now()
while n != 1:
d, stats = (pollard_rho_v0, pollard_rho_v1)[ver](n)
print(d, stats)
assert d > 1, (d, n)
divs.append(d)
assert n % d == 0, (d, n)
n = n // d
pows += min(1, stats['pow'])
print('All divisors:\n', ' '.join(map(str, divs)), sep = '')
print('Avg pow', round(pows / len(divs), 3), ', total time', datetime.datetime.now() - tt)
return divs
p = 1
for i in range(256):
p *= random.randrange(2, 1 << 32)
factor(p, ver = 1)
输出：
................
267890969 {'i': 25551, 'j': 0, 'n_bits': 245, 'd_bits': 28.0, 'pow': 0.523,
'fermat_prp': None, 'time': '0:00:02.363004'}
548977049 {'i': 62089, 'j': 0, 'n_bits': 217, 'd_bits': 29.03, 'pow': 0.5484,
'fermat_prp': None, 'time': '0:00:04.912002'}
3565192801 {'i': 26637, 'j': 0, 'n_bits': 188, 'd_bits': 31.73, 'pow': 0.4633,
'fermat_prp': None, 'time': '0:00:02.011999'}
1044630971 {'i': 114866, 'j': 0, 'n_bits': 156, 'd_bits': 29.96, 'pow': 0.5611,
'fermat_prp': None, 'time': '0:00:06.666996'}
3943786421 {'i': 60186, 'j': 0, 'n_bits': 126, 'd_bits': 31.88, 'pow': 0.4981,
'fermat_prp': None, 'time': '0:00:01.594000'}
3485918759 {'i': 101494, 'j': 0, 'n_bits': 94, 'd_bits': 31.7, 'pow': 0.5247,
'fermat_prp': None, 'time': '0:00:02.161004'}
1772239433 {'i': 102262, 'j': 0, 'n_bits': 63, 'd_bits': 30.72, 'pow': 0.5417,
'fermat_prp': None, 'time': '0:00:01.802996'}
2706462217 {'i': 0, 'j': 1, 'n_bits': 32, 'd_bits': 31.33, 'pow': 0.0,
'fermat_prp': True, 'time': '0:00:00.925801'}
All divisors:
258498 4 99792 121 245864 25 81 2 238008 70 39767 23358624 79 153 27 65 1566 2 31 13 57 1776 446 20 2 3409311 814 37 595384977 2 24 5 147 3738 4514 8372 7 38 237996 430 43 240 1183 10404 11 10234 30 2615625 1263 44590 240 3 101 231 2 79488 799236 2 88059 1578 432500 134 20956 101 3589 155 2471 91 7 6 100608 1995 33 9 181 48 5033 20 16 15 305 44 927 49 76 13 1577 46 144 292 65 2 111890 300 368 430705 6 368 381 1578812 4290 10 48 565 2 2 23606 23020 267 4186 5835 33 4 899 6288 3534 129064 34 315 36190 16900 6 60291 2 12 111631 463 2500 1405 1959 22 112 2 228 3 2192 2 28 321618 4 44 125924200164 9 17956 4224 2848 16 7 162 4 573 843 48 101 224460324 4 768 3 2 8 154 256 2 3 51 784 34 48 14 369 218 9 12 27 152 2 256 2 51 9 9411903 2 131 9 71 6 3 13307904 85608 35982 121669 93 3 3 121 7967 11 20851 19 289 4237 3481 289 89 11 11 121 841 5839 2071 59 29 17293 9367 110801 196219 2136917 631 101 3481 323 101 19 32129 29 19321 19 19 29 19 6113 509 193 1801 347 71 83 1373 191 239 109 1039 2389 1867 349 353 1566871 349 561971 199 1429 373 1231 103 1048871 83 1681 1481 3673 491 691 1709 103 49043 911 673 1427 4147027 569 292681 2153 6709 821 641 569 461 239 2111 2539 6163 3643 5881 2143 7229 593 4391 1531 937 1721 1873 3761 1229 919 178207 54637831 8317 17903 3631 6841 2131 4157 3467 2393 7151 56737 1307 10663 701 2522350423 4253 1303 13009 7457 271549 12391 36131 4943 6899 27077 4943 7723 4567 26959 9029 2063 6607 4721 14563 8783 38803 1889 1613 20479 16231 1847 41131 52201 37507 224351 13757 36299 3457 21739 107713 51169 17981 29173 2287 16253 386611 132137 9181 29123 740533 114769 2287 61553 21121 10501 47269 59077 224951 377809 499729 6257 5903 59999 126823 85199 29501 34589 518113 39409 411667 146603 1044091 312979 291569 158303 41777 115133 508033 154799 13184621 167521 3037 317711 206827 1254059 455381 152639 95531 1231201 494381 237689 163327 651331 351053 152311 103669 245683 1702901 46337 151339 6762257 57787 38959 366343 609179 219749 2058253 634031 263597 540517 1049051 710527 2343527 280967 485647 1107497 822763 862031 583139 482837 1586621 782107 371143 763549 10740361 1372963 62589077 1531627 31991 1206173 678901 4759373 5877959 178439 1736369 687083 53508439 99523 10456609 942943 2196619 376081 802453 10254457 2791597 3231757 2464793 66598351 1535867 16338167 1138639 882953 1483693 12624373 35717041 6427979 5653181 6421873 1434131 1258889 108462803 859667 64298779 261810191 5743483 32314969 5080721 8961767 68011043 7528799 2086957 41618389 19999663 118428929 45556487 40462109 22478363 29039737 17366957 77805557 12775951 50890837 22666991 14892133 691979 133920733 115526921 29092501 2332124099 16835209 101301479 29987047 160734341 35904857 12376361 17774983 2397907 525367681 245240591 48159641 45590383 87274531 69160309 256092673 7430783 588029137 91286513 75817271 393556847 1183839551 71513537 593809903 200299807 161799857 537099259 21510427 335791301 382965337 156133297 180373937 75136921 364790017 174932509 117559207 601612421 54539711 2107325149 566372699 102467207 321156893 1024847609 1250224901 1038888437 3029169139 345512147 4127597891 1043830063 267890969 548977049 3565192801 1044630971 3943786421 3485918759 1772239433 2706462217
Avg pow 0.238 , total time 0:03:48.193658
PS。还决定实施 Pollard-Rho 分解算法的简约但快速版本，纯 Pythonic，准备好复制粘贴到任何项目中（例如分解 Pi 的前 25 位数字）：
def factor(n):
import itertools, math
if n <= 1:
return []
x = 2
for cycle in itertools.count(1):
y = x
for i in range(1 << cycle):
x = (x * x + 1) % n
d = math.gcd(x - y, n)
if d > 1:
return [d] + factor(n // d)
print(factor(1415926535897932384626433))
# [7, 223, 180473, 10739, 468017117899]
我正在尝试使用 python 计算 Pollard rho 数字以获得非常长的整数，例如小于一
65779646778470582047547160396995720887221575959770627441205850493179860146690755880473849736210807898494458426111244201404810495587574110361900128405354081638434164434968839614760264675889940272767106444249
我已经尝试计算我的 intel core i9 10980HK CPU，这导致几分钟的高负载工作没有任何成功。我正在尝试使用 numba 和 @njit 装饰器来连接 RTX 2070 super（在笔记本电脑上），但它给出了以下错误。
- argument 0: Int value is too large:
这里是代码：
import numpy as np
import datetime
def pgcd(a,b):
if b==0:
return a
else:
r=a%b
return pgcd(b,r)
def pollardrho(n):
f = lambda z: z*z+1
x, y, d = 1, 1, 1
c = 0
while d==1:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = pgcd(y-x, n)
return d, c
def test_time(n):
t = datetime.datetime.now()
d, c = pollardrho(int(n))
tps = datetime.datetime.now() - t
print(tps, c, d)
file = open("some_powersmooths_large.txt", "r")
for line in file.readlines():
if not line.startswith("#"):
print(line)
print("\n")
test_time(line)
print("\n")
我该如何处理这种大数计算。
鉴于pollardrho中的操作效率很低，操作需要一段时间我并不奇怪。
但是，我不知道具体的功能，所以我不知道是否可以提高效率。
在Python中，整数具有任意长度。
这意味着它们可以是任意长度，并且 Python 本身将处理使用 64 位整数正确存储它（通过将它们分布在多个整数上）。
（你可以自己测试这个，例如创建一个不能存储在 64 位无符号整数中的整数，比如 a = 2**64，然后检查 a.bit_length() 方法的输出，它应该是 65)
因此，从理论上讲，您应该能够计算任何整数。 但是，由于您使用的是 Numba，因此由于 Numba 的工作方式，您只能使用实际上可以存储在 64 位无符号整数中的整数。 您得到的错误只是数字变得太大而无法存储在 64 位无符号整数中。
底线：没有 Numba，您可以很好地计算出该数字。使用 Numba，你不能。 当然，如果你只想粗略地知道数字是多少，而不是精确地知道，你可以直接使用浮点数。
第 1 部分（共 2 部分，请参阅下面的第 2 部分）。
Numba 最多只能处理 64 位整数，它没有大整数运算，只有 Python 有。 Numba 的开发人员承诺，未来版本将支持大整数。您需要大整数运算，因为您的输入和计算中有非常大的整数。
您的一个优化建议是使用 GMPY2 Python library. It is highly-optimized library of long arithmetics, considerably faster than regular Python implementation of long arithmetics. For very large integers for example it implements multiplication using Fast Fourier Transform，这是目前最快的乘法算法。
但是 GMPY2 安装起来有点困难。 Windows 的最新预编译版本可用 by this link。下载适用于您的 Python 版本的 .whl 文件并通过 pip 安装它，例如对于我的 Windows 64 位 Python 3.7，我下载并安装了 pip install gmpy2-2.0.8-cp37-cp37m-win_amd64.whl。对于 Linux，通过 sudo apt install -y python3-gmpy2.
使用 GMPY2 后，您的代码将变得尽可能快，因为这个库代码几乎是世界上最快的。即使是 Numba（如果它有长算术）也不会改进更多。只有更快的公式和更好的算法才能帮助进一步改进或更小的输入整数。
但是即使使用 GMPY2，您的示例大整数对于您的算法来说也是一种变大的方法。您必须选择更小的整数或更快的算法。我已经 运行 你的算法和数字 5 分钟或更长时间，但没有得到结果。但是，如果使用常规 Python 之前的结果会在 1 小时内完成，那么在使用 GMPY2 之后，它可能会在 10 分钟或更快的时间内完成。
也不是很确定，但可能在你的算法中 f(f(y)) % n 应该等同于 f(f(y) % n) % n 应该计算得更快，因为它会做两次较短的乘法。但这需要额外检查。
此外，您的大整数似乎是质数，正如 Primo elliptic curve based primality proving program, it proved primality of this integer in 3 seconds on my PC. Primo only proves primality (with 100% guarantee) but doesn't factor the number (splitting into divisors). Factoring numbers can be done by programs from this list 所证明的那样，这些程序实现了已知最快的因式分解算法，如果某些链接失效，则 Google 那些程序的名称。
只需将所有整数 n 包装成 gmpy2.mpz(n)。例如，我稍微改进了您的代码，将其包装到 gmpy2.mpz() 中，还制作了一个循环，以便打印所有除数。另外作为一个例子，我没有使用你的大质数，而是使用了一个更小的 - Pi 的前 25 位数字，它是复合的，它的所有除数都在我的 PC 上在 7 秒内打印：
import datetime, numpy as np, gmpy2
def num(n):
return gmpy2.mpz(n)
zero, one = num(0), num(1)
def pgcd(a, b):
if b == zero:
return a
else:
r = a % b
return pgcd(b, r)
def pollardrho(n):
f = lambda z: z * z + one
x, y, d = one, one, one
c = 0
while d == one:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = pgcd(y - x, n)
return d, c
def test_time(n):
n = num(int(n))
divs = []
while n > 1:
t = datetime.datetime.now()
d, c = pollardrho(num(int(n)))
tps = datetime.datetime.now() - t
print(tps, c, d, flush = True)
divs.append(d)
assert n % d == 0, (n, d)
n = n // d
print('All divisors:\n', ' '.join(map(str, divs)), sep = '')
test_time(1415926535897932384626433)
#test_time(65779646778470582047547160396995720887221575959770627441205850493179860146690755880473849736210807898494458426111244201404810495587574110361900128405354081638434164434968839614760264675889940272767106444249)
输出：
0:00:00 2 7
0:00:00 10 223
0:00:00.000994 65 10739
0:00:00.001999 132 180473
0:00:07.278999 579682 468017117899
All divisors:
7 223 10739 180473 468017117899
第 2 部分
阅读维基百科文章 (here and here)，我决定实施更快版本的 Pollard-Rho 算法。
我在下面实现的版本看起来更复杂，但除法和乘法少了两倍，平均而言，循环的总迭代次数也少了。
此改进导致我的测试用例 运行ning 时间 3 分钟，与原始 OP 算法 运行ning 时间 7 分钟相比，我的笔记本电脑。
我的算法也是随机的，这意味着它采用 [1, N - 2] 范围内的随机见证，而不是见证 1 或 2。在极少数情况下，它可能会像维基百科中所说的那样失败，然后我用不同的证人重新 运行 算法。它还使用 Fermat primality test 检查输入数字是否为质数，然后不再搜索任何除数。
为了测试，我使用了由代码 p = 1; for i in range(256): p *= random.randrange(2, 1 << 32) 生成的输入数字 p，基本上它由 256 个因素组成，每个因素最多 32-bits。
我还改进了这两种算法以输出更多统计信息。统计参数之一是 pow，它显示了每个步骤的复杂性，0.25 的 pow 表示如果除数是 d，则当前因式分解步骤花费 c = d^0.25 次迭代来找到该除数 d.正如维基百科中所述，Pollard-Rho 算法平均应具有 pow = 0.25，这意味着找到任何除数 d 的复杂度（迭代次数）约为 d^0.25。
在接下来的代码中还有其他改进，例如提供大量统计信息。并找到循环中的所有因素。
我的测试用例算法版本的平均 pow 为 0.24，原来的先前版本有 0.3。较小的 pow 意味着平均进行较少的循环迭代。
还测试了使用和不使用 GMPY2 的我的版本。显然 GMPY2 对常规 Python 大整数算法没有太大改进，主要是因为 GMPY2 针对真正的大数字（数万位）（使用快速傅立叶变换乘法等）进行了更多优化，而这里的数字并不大在我的测试中。但是 GMPY2 仍然提供了大约 1.35x 倍的加速，提供了 运行 宁时间 3 分钟，而没有 GMPY2 的相同算法几乎 4 分钟。要使用或不使用 gmpy2 进行测试，您只需将 def num(n) 函数内部更改为 return gmpy2.mpz(n) 或 return n.
import datetime, numpy as np, gmpy2, random, math
random.seed(0)
def num(n):
return gmpy2.mpz(n)
zero, one = num(0), num(1)
def gcd(a, b):
while b != zero:
a, b = b, a % b
return a
def pollard_rho_v0(n):
f = lambda z: z * z + one
n, x, y, d, c, t = num(n), one, one, one, 0, datetime.datetime.now()
while d == one:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = gcd(y - x, n)
return d, {'i': c, 'n_bits': n.bit_length(), 'd_bits': round(math.log(d) / math.log(2), 2),
'pow': round(math.log(max(c, 1)) / math.log(d), 4), 'time': str(datetime.datetime.now() - t)}
def is_fermat_prp(n, trials = 32):
n = num(n)
for i in range(trials):
a = num((3, 5, 7)[i] if i < 3 else random.randint(2, n - 2))
if pow(a, n - 1, n) != 1:
return False
return True
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
# https://ru.wikipedia.org/wiki/%D0%A0%D0%BE-%D0%B0%D0%BB%D0%B3%D0%BE%D1%80%D0%B8%D1%82%D0%BC_%D0%9F%D0%BE%D0%BB%D0%BB%D0%B0%D1%80%D0%B4%D0%B0
def pollard_rho_v1(N):
AbsD = lambda a, b: a - b if a >= b else b - a
N, fermat_prp, t = num(N), None, datetime.datetime.now()
SecsPassed = lambda: (datetime.datetime.now() - t).total_seconds()
for j in range(8):
i, stage, y, x = 0, 2, num(1), num(random.randint(1, N - 2))
while True:
if (i & 0x3FF) == 0 and fermat_prp is None and (SecsPassed() >= 15 or j > 0):
fermat_prp = is_fermat_prp(N)
if fermat_prp:
r = N
break
r = gcd(N, AbsD(x, y))
if r != one:
break
if i == stage:
y = x
stage <<= one
x = (x * x + one) % N
i += 1
if r != N or fermat_prp:
return r, {'i': i, 'j': j, 'n_bits': N.bit_length(), 'd_bits': round(math.log(r) / math.log(2), 2),
'pow': round(math.log(max(i, 1)) / math.log(r), 4), 'fermat_prp': fermat_prp, 'time': str(datetime.datetime.now() - t)}
assert False, f'Pollard-Rho failed after {j + 1} trials! N = {N}'
def factor(n, *, ver = 1):
assert n > 0, n
n, divs, pows, tt = int(n), [], 0., datetime.datetime.now()
while n != 1:
d, stats = (pollard_rho_v0, pollard_rho_v1)[ver](n)
print(d, stats)
assert d > 1, (d, n)
divs.append(d)
assert n % d == 0, (d, n)
n = n // d
pows += min(1, stats['pow'])
print('All divisors:\n', ' '.join(map(str, divs)), sep = '')
print('Avg pow', round(pows / len(divs), 3), ', total time', datetime.datetime.now() - tt)
return divs
p = 1
for i in range(256):
p *= random.randrange(2, 1 << 32)
factor(p, ver = 1)
输出：
................
267890969 {'i': 25551, 'j': 0, 'n_bits': 245, 'd_bits': 28.0, 'pow': 0.523,
'fermat_prp': None, 'time': '0:00:02.363004'}
548977049 {'i': 62089, 'j': 0, 'n_bits': 217, 'd_bits': 29.03, 'pow': 0.5484,
'fermat_prp': None, 'time': '0:00:04.912002'}
3565192801 {'i': 26637, 'j': 0, 'n_bits': 188, 'd_bits': 31.73, 'pow': 0.4633,
'fermat_prp': None, 'time': '0:00:02.011999'}
1044630971 {'i': 114866, 'j': 0, 'n_bits': 156, 'd_bits': 29.96, 'pow': 0.5611,
'fermat_prp': None, 'time': '0:00:06.666996'}
3943786421 {'i': 60186, 'j': 0, 'n_bits': 126, 'd_bits': 31.88, 'pow': 0.4981,
'fermat_prp': None, 'time': '0:00:01.594000'}
3485918759 {'i': 101494, 'j': 0, 'n_bits': 94, 'd_bits': 31.7, 'pow': 0.5247,
'fermat_prp': None, 'time': '0:00:02.161004'}
1772239433 {'i': 102262, 'j': 0, 'n_bits': 63, 'd_bits': 30.72, 'pow': 0.5417,
'fermat_prp': None, 'time': '0:00:01.802996'}
2706462217 {'i': 0, 'j': 1, 'n_bits': 32, 'd_bits': 31.33, 'pow': 0.0,
'fermat_prp': True, 'time': '0:00:00.925801'}
All divisors:
258498 4 99792 121 245864 25 81 2 238008 70 39767 23358624 79 153 27 65 1566 2 31 13 57 1776 446 20 2 3409311 814 37 595384977 2 24 5 147 3738 4514 8372 7 38 237996 430 43 240 1183 10404 11 10234 30 2615625 1263 44590 240 3 101 231 2 79488 799236 2 88059 1578 432500 134 20956 101 3589 155 2471 91 7 6 100608 1995 33 9 181 48 5033 20 16 15 305 44 927 49 76 13 1577 46 144 292 65 2 111890 300 368 430705 6 368 381 1578812 4290 10 48 565 2 2 23606 23020 267 4186 5835 33 4 899 6288 3534 129064 34 315 36190 16900 6 60291 2 12 111631 463 2500 1405 1959 22 112 2 228 3 2192 2 28 321618 4 44 125924200164 9 17956 4224 2848 16 7 162 4 573 843 48 101 224460324 4 768 3 2 8 154 256 2 3 51 784 34 48 14 369 218 9 12 27 152 2 256 2 51 9 9411903 2 131 9 71 6 3 13307904 85608 35982 121669 93 3 3 121 7967 11 20851 19 289 4237 3481 289 89 11 11 121 841 5839 2071 59 29 17293 9367 110801 196219 2136917 631 101 3481 323 101 19 32129 29 19321 19 19 29 19 6113 509 193 1801 347 71 83 1373 191 239 109 1039 2389 1867 349 353 1566871 349 561971 199 1429 373 1231 103 1048871 83 1681 1481 3673 491 691 1709 103 49043 911 673 1427 4147027 569 292681 2153 6709 821 641 569 461 239 2111 2539 6163 3643 5881 2143 7229 593 4391 1531 937 1721 1873 3761 1229 919 178207 54637831 8317 17903 3631 6841 2131 4157 3467 2393 7151 56737 1307 10663 701 2522350423 4253 1303 13009 7457 271549 12391 36131 4943 6899 27077 4943 7723 4567 26959 9029 2063 6607 4721 14563 8783 38803 1889 1613 20479 16231 1847 41131 52201 37507 224351 13757 36299 3457 21739 107713 51169 17981 29173 2287 16253 386611 132137 9181 29123 740533 114769 2287 61553 21121 10501 47269 59077 224951 377809 499729 6257 5903 59999 126823 85199 29501 34589 518113 39409 411667 146603 1044091 312979 291569 158303 41777 115133 508033 154799 13184621 167521 3037 317711 206827 1254059 455381 152639 95531 1231201 494381 237689 163327 651331 351053 152311 103669 245683 1702901 46337 151339 6762257 57787 38959 366343 609179 219749 2058253 634031 263597 540517 1049051 710527 2343527 280967 485647 1107497 822763 862031 583139 482837 1586621 782107 371143 763549 10740361 1372963 62589077 1531627 31991 1206173 678901 4759373 5877959 178439 1736369 687083 53508439 99523 10456609 942943 2196619 376081 802453 10254457 2791597 3231757 2464793 66598351 1535867 16338167 1138639 882953 1483693 12624373 35717041 6427979 5653181 6421873 1434131 1258889 108462803 859667 64298779 261810191 5743483 32314969 5080721 8961767 68011043 7528799 2086957 41618389 19999663 118428929 45556487 40462109 22478363 29039737 17366957 77805557 12775951 50890837 22666991 14892133 691979 133920733 115526921 29092501 2332124099 16835209 101301479 29987047 160734341 35904857 12376361 17774983 2397907 525367681 245240591 48159641 45590383 87274531 69160309 256092673 7430783 588029137 91286513 75817271 393556847 1183839551 71513537 593809903 200299807 161799857 537099259 21510427 335791301 382965337 156133297 180373937 75136921 364790017 174932509 117559207 601612421 54539711 2107325149 566372699 102467207 321156893 1024847609 1250224901 1038888437 3029169139 345512147 4127597891 1043830063 267890969 548977049 3565192801 1044630971 3943786421 3485918759 1772239433 2706462217
Avg pow 0.238 , total time 0:03:48.193658
PS。还决定实施 Pollard-Rho 分解算法的简约但快速版本，纯 Pythonic，准备好复制粘贴到任何项目中（例如分解 Pi 的前 25 位数字）：
def factor(n):
import itertools, math
if n <= 1:
return []
x = 2
for cycle in itertools.count(1):
y = x
for i in range(1 << cycle):
x = (x * x + 1) % n
d = math.gcd(x - y, n)
if d > 1:
return [d] + factor(n // d)
print(factor(1415926535897932384626433))
# [7, 223, 180473, 10739, 468017117899]