在Python中使用gmpy计算大数时崩溃？

Question

我被推荐使用gmpy来帮助高效计算大数字。在此之前，我只是用Python写代码，结果我的程序运行了一两天后就内存不够用了（我也不太明白为什么会这样，因为我的程序内存使用应该是基本保持不变的……可能是内存泄漏？）

总之，我在运行程序几秒钟后就一直遇到一个奇怪的错误：

mp_allocate< 545275904->545275904 >
Fatal Python error: mp_allocate failure

This application has requested the Runtime to terminate it in an unusual way. 
Please contact the application's support team for more information.

而且，Python崩溃了，Windows 7给我弹出了一个通用的提示框，显示python.exe已停止工作。

使用标准的Python整数时没有发生这种情况。现在我切换到gmpy后，运行脚本几秒钟就出现这个错误。我以为gmpy是专门用来处理大数字运算的呢？

为了参考，这里有一个会产生错误的示例程序：

import gmpy2

p = gmpy2.xmpz(3000000000)
s = gmpy2.xmpz(2)
M = s**p

for x in range(p):
    s = (s * s) % M

我有10GB的内存，使用gmpy之前这个脚本可以运行好几天而不会内存不够（我还是不太明白为什么会这样，因为s的大小其实并没有变得很大……）

有没有人有什么想法？

编辑：我忘了提，我使用的是Python 3.2

Answer 1

免责声明：我是 gmpy 和 gmpy2 的维护者。

我今晚才能测试这个。不过我有几个评论和问题。

与其用 (s * s) % M，不如用 pow(s, 2, M)。这样应该会更快。

如果你用 gmpy2.mpz() 而不是 gmpy2.xmpz() 会发生什么？

你是在运行 64 位版本的 Python 和 gmpy2 吗？（我猜是这样，但我想确认一下。）

关于 range 和 xrange，在 Python 3.x 中，range 替代了 xrange。

编辑补充信息。

崩溃的原因是因为在 32 位版本中内部结构溢出。使用 64 位版本的 Python 和 gmpy 或 gmpy2 是正确的解决办法。

unpack(x,n) 函数类似于字符串的 split()：它把一个数字分成一系列 n 位的值。它的效果和下面的代码相同，但速度更快：

def unpack(x,n):
r = []
m = 2**n
while x:
    x, temp = divmod(x,m)
    r.append(temp)
return r

一些文档可以通过 help(gmpy2.unpack) 获取，但更好的文档在我的待办事项列表上。

unpack() 能够用来消除 % 操作的原因，和检查一个十进制数字是否能被 9 整除时加各位数字的原理是一样的。在这个例子中，unpack() 创建了 p 位的数字，而我们是用 2**p - 1 来进行除法。

这里有一些测试代码：

import gmpy2
import time

def mersenne1(p):
    '''Primality test for Mersenne prime: 2**p -1.
    Uses native Python longs. Does not verify that p is prime.'''

    s = 4
    M = 2**p - 1
    for i in range(p-2):
        s = ((s*s)-2) % M
    return False if s else True

def mersenne2(p):
    '''Primality test for Mersenne prime: 2**p -1.
    Uses gmpy2.mpz. Does not verify that p is prime.'''

    s = gmpy2.mpz(4)
    M = gmpy2.mpz(2)**p - 1
    for i in range(p-2):
        s = ((s*s)-2) % M
    return False if s else True

def mersenne3(p):
    '''Primality test for Mersenne prime: 2**p -1.
    Uses gmpy2.mpz and no mod. Does not verify that p is prime.'''

    s = gmpy2.mpz(4)
    M = gmpy2.mpz(2)**p - 1
    for i in range(p-2):
        s = (s*s)
        s = sum(gmpy2.unpack(s, p))
        s = sum(gmpy2.unpack(s, p))
        if s < 2:
            s = M - 2 + s
        else:
            s = s - 2
    return False if s else True

if __name__ == "__main__":
    p = 44497

    start = time.time()
    result1 = mersenne1(p)
    print("Elapsed time: {:6.3f}".format(time.time() - start))

    start = time.time()
    result2 = mersenne2(p)
    print("Elapsed time: {:6.3f}".format(time.time() - start))

    start = time.time()
    result3 = mersenne3(p)
    print("Elapsed time: {:6.3f}".format(time.time() - start))

    if result1 == result2 == result3:
        print("All three tests are equal!")
    else:
        print("Oops, something has gone wrong.")

还有一些运行时间...

C:\x64\Python32>python.exe mersenne.py
Elapsed time: 163.683
Elapsed time: 12.782
Elapsed time:  3.630
All three tests are equal!