问题描述
我正在尝试使用 python 计算 Pollard rho 数，用于非常长的整数，例如低于 1
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() 方法的输出，它应该说 { {1}})
所以，从理论上讲，您应该能够计算任何整数。 但是，由于您使用的是 Numba，因此由于 Numba 的工作方式，您只能使用可以实际存储在 64 位无符号整数中的整数。 您得到的错误只是数字变得太大而无法存储在 64 位无符号整数中。
底线：如果没有 Numba，您可以很好地计算该数字。使用 Numba，您不能。 当然，如果你只是想知道大概的数字是多少，而不是精确的，你可以用浮点数代替。
,第 1 部分（共 2 个，请参阅下面的第 2 部分）。
Numba 最多只能处理 64 位整数，它没有大整数算法，只有 Python 有。由于 Numba 的开发人员承诺，未来版本将支持大整数。您需要大整数算术，因为您的输入和计算中有非常大的整数。
给您的一个优化建议是使用 GMPY2 Python 库。它是高度优化的长算术库，比长算术的常规 Python 实现快得多。例如，对于非常大的整数，它使用 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 椭圆曲线素数证明程序所证明的那样，它在我的 PC 上在 3 秒内证明了这个整数的素数。 Primo 只证明素数（100% 保证），但不分解数字（分成除数）。因数可以由程序 from this list 完成，这些程序实现了最快的已知因数算法，如果某些链接失效，则谷歌这些程序名称。
只需将所有整数 n 包装成 gmpy2.mpz(n)。例如，我稍微改进了您的代码，将其包装到 gmpy2.mpz() 中并进行了循环，以便打印所有除数。同样作为一个例子，我没有取你的大素数，而是取一个小得多的 - Pi 的前 25 位数字，它是复合的，它的所有除数都在 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 和 here）后，我决定实施一个更快版本的 Pollard-Rho 算法。
我在下面实现的版本看起来更复杂，但除法和乘法少了两倍，而且平均而言，循环的总迭代次数也更少。
与在笔记本电脑上运行时间为 3 分钟的原始 OP 算法相比，这种改进导致我的测试用例的运行时间为 7 分钟。
我的算法也是随机的，这意味着它采用 [1,N - 2] 范围内的随机见证而不是见证 1 或 2。在极少数情况下，它可能会像维基百科中所说的那样失败，然后我用不同的见证人重新运行算法。它还使用 Fermat primality test 来检查输入数是否为质数，然后不再搜索任何除数。
对于测试，我使用由代码 p 生成的输入编号 p = 1; for i in range(256): p *= random.randrange(2,1 << 32)，基本上它由 256 个因子组成，每个 32-bits 最多。
我还改进了两种算法以输出更多统计信息。其中一个统计参数是 pow，它显示了每个步骤的复杂性，0.25 的 pow 表明如果除数是 d，那么当前的分解步骤花费了 c = d^0.25 次迭代来找到这个除数 {{1} }.正如维基百科中所说的 Pollard-Rho 算法平均应该有 d，这意味着找到任何除数 pow = 0.25 的复杂度（迭代次数）大约是 d。
在下一个代码中，还有其他改进，例如在途中提供了大量统计信息。并找到循环中的所有因素。
我的测试用例的算法版本的平均 pow 为 d^0.25，原始以前的版本为 0.24。较小的 pow 意味着平均执行较少的循环迭代。
还测试了我的版本，有和没有 GMPY2。显然 GMPY2 对常规 Python 大整数算法的改进不大，主要是因为 GMPY2 对真正的大数字（数万位）（使用快速傅立叶变换乘法等）进行了更优化，而在我的测试中，这里的数字并不大。但是 GMPY2 仍然提供了大约 0.3 倍的加速，提供了 1.35x 分钟的运行时间，而没有 GMPY2 的相同算法的运行时间几乎为 3 分钟。要使用或不使用 gmpy2 进行测试，您只需将 4 函数内部更改为 def num(n) 或 return gmpy2.mpz(n)。
return n输出：
import datetime,gmpy2,random,math
random.seed(0)
def num(n):
return gmpy2.mpz(n)
zero,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,t = num(n),datetime.datetime.now()
while d == one:
c += 1
x = f(x) % n
y = f(f(y)) % n
d = gcd(y - x,{'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,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),'pow': round(math.log(max(i,1)) / math.log(r),'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,n)
n = n // d
pows += min(1,stats['pow'])
print('All divisors:\n',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)
附注。还决定实现简约但快速的 Pollard-Rho 分解算法版本，纯 Pythonic，准备复制粘贴到任何项目中（例如分解 Pi 的前 25 位数字）：
................
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,'n_bits': 217,'d_bits': 29.03,'pow': 0.5484,'time': '0:00:04.912002'}
3565192801 {'i': 26637,'n_bits': 188,'d_bits': 31.73,'pow': 0.4633,'time': '0:00:02.011999'}
1044630971 {'i': 114866,'n_bits': 156,'d_bits': 29.96,'pow': 0.5611,'time': '0:00:06.666996'}
3943786421 {'i': 60186,'n_bits': 126,'d_bits': 31.88,'pow': 0.4981,'time': '0:00:01.594000'}
3485918759 {'i': 101494,'n_bits': 94,'d_bits': 31.7,'pow': 0.5247,'time': '0:00:02.161004'}
1772239433 {'i': 102262,'n_bits': 63,'d_bits': 30.72,'pow': 0.5417,'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
设置时间 控制面板
错误1：Request method ‘DELETE‘ not supported 错误还原：...