在Python中使用GMP（gmpy2）

之前笔者曾写过《初试在Python中使用PARI/GP》，简单介绍了一下在Python中调用PARI/GP的方法。PARI/GP是一个比较强大的数论库，“针对数论中的快速计算（大数分解，代数数论，椭圆曲线...）而设计”，它既可以被C/C++或Python之类的编程语言调用，而且它本身又是一种自成一体的脚本语言。而如果仅仅需要高精度的大数运算功能，那么GMP似乎更满足我们的需求。

了解C/C++的读者都会知道GMP（全称是GNU Multiple Precision Arithmetic Library，即GNU高精度算术运算库），它是一个开源的高精度运算库，其中不但有普通的整数、实数、浮点数的高精度运算，还有随机数生成，尤其是提供了非常完备的数论中的运算接口，比如Miller-Rabin素数测试算法、大素数生成、欧几里德算法、求域中元素的逆、Jacobi符号、legendre符号等^[来源]。虽然在C/C++中调用GMP并不算复杂，但是如果能在以高开发效率著称的Python中使用GMP，那么无疑是一件快事。这正是本文要说的gmpy2。

gmpy2介绍

gmpy2是Python的一个扩展库，是对GMP的封装，它的前身是gmpy。经过作者的调整和封装，使得gmpy2的使用大大简化。根据笔者目前的体会，在Python中使用gmpy2，跟在C/C++中使用GMP相比，至少以下优势：

1、简化的函数名：作者对GMP的大量函数名做了简化，比如素数的概率性测试，在C/C++中，命令是mpz_probab_prime_p，在gmpy2中，只需要简单的is_prime；

2、方便的运算符重载：如果在C/C++中，两个mpz整数相加，我们要用mpz_add(z_i, z_i, z_o)，在gmpy2中，我们只需要用+号即可，更一般的，将一个mpz整数与int整数相加，也只需要+号，类似的还有减、乘、除、求模、乘方等。

笔者所学还不深入，因此只能这样片面的评论了，当然，作为第三方语言的调用，gmpy2的效率总比直接在C/C++中直接调用GMP会低一些，但差别是很小的，因为gmpy2就是预先编译好的C语言库而已。具体的gmpy2教程请看：
http://gmpy2.readthedocs.org/en/latest/

基本使用

本文只做简单介绍。以下代码均在Python 3.4中运行。初始化一个大整数，只需要

import gmpy2
n=gmpy2.mpz(1257787) #初始化
gmpy2.is_prime(n) #概率性素性测试

这里跟C/C++是平行的，其实括号里边的参数，可以是整型，也可以是字符串。gmpy2中不仅集成了大整数mpz，还有高精度浮点数mpfr，用法跟C/C++都是平行的。下面是一些基本计算

a+2 #求和，结果是一个mpz
a-2 #求差，结果是一个mpz
a*2 #求积，结果是一个mpz
a/2 #求商，结果是一个mpfr
a//2 #求商，结果是一个mpz
a**2 #求平方，结果是一个mpz
a%2 #求模，结果是一个mpz

这些运算符的重载既直观又方便，其中将数字2换成mpz类型变量也是可以的，这自然不用多说。事实上，运算a+2，就是先将2转换为mpz中的2，然后执行mpz中的加法。

与Python自身比较

Python本身支持高精度大整数运算，一般的数字比较小的范围足够了，只是并非“快速”，读者只需要分别运行以下两句代码

gmpy2.mpz(1257787)**123456 #gmpy2计算
1257787**123456 #Python自带计算

就能感觉到了（当然，可能读者的配置还不错，两者并没有明显差别，那么请把指数扩大十倍~笔者的笔记本配置不高，见笑了^_^）。

尝试大数分解

下面来一段自己编写的大数分解的代码：

from gmpy2 import *
import time
start=time.clock()
n = mpz(63281217910257742583918406571)
x = mpz(2)
y = x**2 + 1
for i in range(n):
    p = gcd(y-x,n)
    if p != 1:
        print(p)
        break
    else:
        y=(((y**2+1)%n)**2+1)%n
        x=(x**2+1)%n
end=time.clock()
print(end-start)

该代码调用了gmpy2，算法则是Pollard rho方法，只是最单纯地运行，并没做任何优化。该算法在我的笔记本上分解

63281217910257742583918406571 = 125778791843321 × 503115167373251

耗时84秒（获得一个因子即停止）。当然这不是什么了不起的成绩，这个数字在Mathematica和PARI/GP中都是瞬间完成的。这个脚本只是纯粹练习和测试Python对gmpy2的调用。

关于Pollard rho方法的步骤：

1、给定合数$n$，以及$x_0=2,y_0=x_0^2+1$；

2、计算$p=Gcd(x_0-y_0,n)$，如果p非1，则结果为$n$的一个因子，停止；

3、如果p=1，则按照$x_{n+1}=x_n^2+1,y_{n+1}=(y_n^2+1)^2+1$，然后重复第2步。

其中，为了减少计算量，每一步$x_n ,y_n$的计算，都要作模$n$运算，以降低计算量。上述步骤只是最简单的计算部分，还没有处理各种可能出现的异常，比如出现$n|(x_n -y_n)$等；也没有对某些细节进行优化，请谨慎引用。

参考网址：
http://hi.baidu.com/pytvzcnbuolrsue/item/ae714592f1fda8d87b7f010a

转载到请包括本文地址：https://kexue.fm/archives/3026

更详细的转载事宜请参考：《科学空间FAQ》