BoJone很喜欢Python,也很喜欢数论,所以就喜欢利用Python玩数论了。平时也喜欢自己动手写一些数论函数,毕竟Python支持大整数高精度运算,这点是非常好的;但是,在很多实际应用中,还是希望能有一个现成的数论函数库来调用。之前尝试过数学研发网的HugeCalc库,但是由于各种不熟悉不了了之。后来论坛上的无心老兄推荐了PARI/GP,小试一下,居然在Python上成功调用了。以后再也不用担心Python上的数论计算问题了,呵呵~

先来一段官方的关于PARI/GP的介绍:

PARI/GP 是一种针对数论中的快速计算(大数分解,代数数论,椭圆曲线...) 而设计的广泛应用的计算机代数系统,同样具备大量实用的函数来对于数学实体的计算, 诸如矩阵,多项式,幂级数,代数数,以及相当多的超越方程等等。 PARI也可以作为快速计算的C语言库。

这个系统最初是由Henri Cohen和他的合作者们(来自法国的波尔多第一大学)发 展起来的。 PARI现在在公共授权条款(GPL)之下并在许多志愿者的帮助下由Karim Belabas进行维护。

-PARI 是一个可以快速计算的C语言库。
-gp是一个容易上手的交互式外壳,它赋予用户调用PARI函数的权力。
-GP 是gp脚本语言的名字。
-gp2c, GP到C的编译器,通过编译GP脚本至C语言以及透明载入结果函数至gp来结合两者的最优点。 (由gp2c编译的脚本一般情况下运行速度将快3至4倍。) 目前gp2c只能运行GP语言的一部分。

PARI/GP提供了源码包,可以编译使用,也提供了Windows下的可执行安装包,直接安装使用也可。本文只讨论在Windows平台下用Python调用PARI,Linux平台估计可以类似地调用,但写这篇文章之前,我还没亲自测试过。Windows下的调用只需要PARI/GP安装目录下的libpari.dll文件,读者可以从官网下载安装包安装后提取它,也可以在本文末尾的附件中下载它。

首先,需要注意的是,调用PARI的Python版本必须是32位的,Python 2和Python 3均可,但64位会出错(可以在64位的操作系统中,安装32位的Python),当然,也可以用PyPy来调用。PyPy基于32位的Python,最新版本已经支持Python 3,读者可以选择自己喜欢的。先看在Python 2.x中的方法,3.x在细节上有些差别,在后面会说明。

DLL是动态链接库,加载DLL之后,就可以使用DLL里边的函数。这跟加载模块有点类似。不同的是,DLL一般是用C语言写的,运行效率比较高,而在Python中加载C/C++的函数,可以发挥C/C++的运行效率和Python的开发效率。在Python中加载DLL,需用导入ctypes模块,导入之后,利用ctypes.cdll.LoadLibrary()函数载入libpari.dll(linux下也有类似的链接库,后缀为.so)。代码如下

import ctypes
from ctypes import * 
pari = ctypes.cdll.LoadLibrary('libpari.dll')

在调用PARI的功能之前,首先得用它的pari_init()函数初始化,比如pari.pari_init(4000000,2),其中第一个参数表示提供给PARI使用的字节数,也就是PARI最大内存占用,一般来说不要小于500000,第二个参数是预先生成的素数表,预先生成部分素数表可以在之后处理素数相关问题时,降低一点计算量,但是这个不是必须的,你也可以设它为0,它同样还是可以完成素数相关运算。

下面以使用PARI的primes()函数为例。prime(n)函数的作用是输出前n个素数。代码如下

import ctypes
from ctypes import * 
pari = ctypes.cdll.LoadLibrary('libpari.dll') #导入libpari.dll

pari.pari_init(4000000,2) #初始化
pari.primes.restype = ctypes.POINTER(ctypes.c_long) #这两句不大清楚作用
pari.GENtostr.restype = ctypes.POINTER(ctypes.c_char) #这两句不大清楚作用

n = pari.primes(10) #生成
y = pari.GENtostr(n) #转换
print(ctypes.string_at(y)) #输出

由于只是刚刚尝试,我也不是很懂restype的作用,望文生义,应该就是把函数的类型跟指针相对应。这代码是参考网络上的一些代码而来的,我测试也发现,就算没有restype的这两句,程序也能够正常工作,不知道有什么不同。程序的过程是先用primes()函数生成结果,结果是GEN类型的,在Python中没法识别,用GENtostr()函数转换为字符串,但这函数返回的是地址,要用string_at()读取输出。在Python 2.7中,运行结果为(输出的是字符串):

[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]

以上代码在Python 3同样可以正常运行,但是输出结果有小许不同。以下会说明。

第二个例子是isprime()函数,这是素性测试函数。这种函数有些特殊,它不是用C/C++的int类型,因此就不能简单输入Python的整数,估计PARI自定义了自己的大整数类型,或者用其他的方式储存大整数。我们要从字符串中读入数字,代码如下:

import ctypes
from ctypes import * 
pari = ctypes.cdll.LoadLibrary('libpari.dll') #导入libpari.dll

pari.pari_init(4000000,2) #初始化
pari.gp_read_str.restype = ctypes.POINTER(ctypes.c_int) #这句不大清楚作用

n=pari.isprime(pari.gp_read_str(str(127))) #测试127的素性
print(n)

核心代码就是倒数第二句,首先将127转化为字符串,然后用gp_read_str()函数从字符串读入,最后用isprime()函数判断,本例的输出结果为1。

但是,如果把这个例子放到Python 3中运行,就会出现错误。原因在于,Python 3的字符串str,不等于Python 2.x的str了,在3中定义了一种新的类型bytes,事实上,2.x的str在这里更像3中的bytes,所以,在Python中,代码要改成

pari.isprime(pari.gp_read_str(b'127'))

或者(用str的encode()函数转成bytes类型)

pari.isprime(pari.gp_read_str(str(127).encode()))

同样,在Python 3中,ctypes.string_at(y)的返回对象也不是str,而是bytes,读者可以发现,输出结果多带了一个b,正是这个原因。

暂时就介绍到这里了,如果有新的发现,会继续跟大家分享~~希望高手路过多多指教。

下载:libpari.zip

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

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

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (Jul. 22, 2014). 《初试在Python中使用PARI/GP 》[Blog post]. Retrieved from https://www.kexue.fm/archives/2775

@online{kexuefm-2775,
        title={初试在Python中使用PARI/GP},
        author={苏剑林},
        year={2014},
        month={Jul},
        url={\url{https://www.kexue.fm/archives/2775}},
}