如何在Python中计算很长的数字

Question

我正在尝试使用python计算Pollard rho数，计算非常长的整数，例如小于1的整数

65779646778470582047547160396995720887221575959770627441205850493179860146690755880473849736210807898494458426111244201404810495587574110361900128405354081638434164434968839614760264675889940272767106444249

我已经尝试在我的intel core i9 10980HKCPU上进行计算，这会导致几分钟的高负载工作，但没有任何成功。我试图使用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")

我如何处理这种类型的大数计算

Answer 1

第1部分（第2部分，见下文第2部分）

Numba最多只能处理64位整数，它没有大整数算法，只有Python有。作为Numba promise的开发者，大整数将在未来的版本中得到支持。您需要大整数算术，因为您的输入和计算中有非常大的整数

一个优化建议是使用GMPY2Python库。它是高度优化的长算术库，比长算术的常规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椭圆曲线的素数证明程序所证明的，它在我的电脑上3秒钟内证明了这个整数的素数。Primo只证明素数（100%保证），但不计算数（分解为除数）。分解数字可以由程序from this list完成，这些程序实现已知最快的分解算法，如果一些链接失效，那么用谷歌搜索这些程序的名称

只需将所有整数n包装成gmpy2.mpz(n)。例如，我对您的代码进行了一些改进，将其包装到gmpy2.mpz()中，并创建了一个循环，以便打印所有除数。同样作为一个例子，我不是拿你的大素数，而是一个小得多的Pi的前25位数字，它是复合的，它的所有除数都在我的电脑上以7秒的时间打印出来：

Try it online!

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

---

第二部分

阅读维基百科的文章（here和here），我决定实现一个更快版本的Pollard Rho算法

下面实现的我的版本看起来更复杂，但除法和乘法减少了两倍，平均而言，循环的迭代次数也减少了一倍

与我的笔记本电脑上运行时间为7分钟的原始OP算法相比，这种改进导致我的测试用例的运行时间为3分钟

我的算法也是随机的，这意味着它不使用见证1或见证2，而是在范围[1, N - 2]内使用随机见证。在极少数情况下，它可能会像维基百科中所说的那样失败，然后我用不同的见证人重新运行算法。此外，它还使用Fermat primality test检查输入的数字是否为素数，然后不搜索更多的除数

对于测试，我使用了由代码p = 1; for i in range(256): p *= random.randrange(2, 1 << 32)生成的输入数p，基本上它最多由256个因子组成

我还改进了这两种算法以输出更多的统计数据。其中一个统计参数是pow，它显示了每个步骤的复杂性，pow为0.25表示如果除数是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

Try it online!

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位数字）：

Try it online!

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]

Answer 2

考虑到pollardrho中的操作非常低效，我对该操作需要一段时间并不感到惊讶。 然而，我不知道这个特定的函数，所以我不知道它是否可以变得更有效

在Python中，整数具有任意长度。 这意味着它们可以是任意长度的，Python本身将使用64位整数（通过将它们分散到多个整数上）正确地存储它们。 （您可以自己测试一下，例如创建一个不能存储在64位无符号整数中的整数，如a = 2**64，然后检查a.bit_length()方法的输出，该方法应该是65）

所以，从理论上讲，你应该能够计算任何整数。 但是，由于您使用的是Numba，因此由于Numba的工作方式，您仅限于可以实际存储在64位无符号整数中的整数。 您得到的错误只是数字变得太大，无法存储在64位无符号整数中

底线：没有Numba，你可以很好地计算这个数字。有了麻木，你就不能。 当然，如果您只想大致知道数字是什么，而不想知道精确的数字，那么您可以使用浮点数