如何使用Tuple/Array/Vector从Python(ctypes)调用PARI/GP? [英] How to use Tuple/Array/Vector to call PARI/GP from Python (ctypes)?

查看:110
本文介绍了如何使用Tuple/Array/Vector从Python(ctypes)调用PARI/GP?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我想从Python调用 PARI/GP .我需要使用PARI的ellisdivisible(E; P; n;{&Q})函数 (请参阅此链接中第441页的函数3.15.35 :) ,所以我必须传递2个向量或数组(例如E = ellinit([0,-1,1,0,0], K);P = [0,0];),我该怎么做?

要从Python (由Thomas Baruchel提供)中调用单参数/变量的PARI函数(用C语言编写),我们有类似下面的内容-

import ctypes

# load the library
pari=ctypes.cdll.LoadLibrary("libpari.so")

# set the right return type of the functions
pari.stoi.restype = ctypes.POINTER(ctypes.c_long)
pari.nextprime.restype = ctypes.POINTER(ctypes.c_long)

# initialize the library 
pari.pari_init(2**19,0)

def nextprime(v):
  g = pari.nextprime(pari.stoi(ctypes.c_long(v))) # nextprime(argument) is a PARI function
  return pari.itos(g)



print( nextprime(456) )

例如,我尝试过-

h=(0,0,0, 4,6)
pari.stoi.restype = ctypes.POINTER(ctypes.c_long*5)
pari.ellinit.restype = ctypes.POINTER(ctypes.c_long)
def ellinit(v):
  g = pari.ellinit(pari.stoi(ctypes.c_long(v)*5))
  return pari.itos(g)


print(ellinit(h))

我遇到了以下错误-

   File "C:\Users\miron\Desktop\trash5\x\f.py", line 68, in <module>
    print( ellinit(h) )
  File "C:\Users\miron\Desktop\trash5\x\f.py", line 62, in ellinit
    g = pari.ellinit(pari.stoi(ctypes.c_long(v)*5))
TypeError: an integer is required (got type tuple)

如何传递元组/数组/向量?谢谢.

修改: 尝试获取ellisdivisible(E; P; n;{&Q})-

失败

from ctypes import *

pari = cdll.LoadLibrary("C:\\Program Files\\Python37\\Pari64-2-11-3\\libpari.dll")

pari.stoi.restype = POINTER(c_long)
pari.cgetg.restype = POINTER(POINTER(c_long))
pari.ellinit.restype = POINTER(POINTER(c_long))
#-------------------------CHANGE 1
pari.ellisdivisible.restype = c_long
Q = pari.stoi(c_long(0))
#-------------------------
(t_VEC, t_COL, t_MAT) = (17, 18, 19)  # incomplete
precision = c_long(38)
pari.pari_init(2 ** 19, 0)
def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1
def main():
    h = (0, 0, 0, 0, 1)
    P=(0,0)
    res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision)
#---------------CHANGE 2
   # res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision).disc
    y = pari.ellisdivisible(res, t_vec(P), pari.stoi(c_long(5)), byref(Q))
    print(pari.itos(y))
#---------------
    for i in range(1, 13):
        print(pari.itos(res[i]))

if __name__ == '__main__':
    main()

错误是-

Traceback (most recent call last):
  File "C:\Users\miron\Desktop\trash5\x\ex - Copy (2).py", line 34, in <module>
    main()
  File "C:\Users\miron\Desktop\trash5\x\ex - Copy (2).py", line 28, in main
    print(pari.itos(y))
OSError: exception: access violation reading 0x0000000000000009

解决方案

Python元组或C数组无法直接使用,因为PARI使用的是PARI/GP特定矢量,其类型/长度在第一个元素中进行了编码.

在第4.4.1节Creation of PARI objects中说:

创建PARI对象的基本功能是GEN cgetg(long l,long t)l以符号形式指定要分配给该对象的长字数,并且t是该对象的类型(请参见第4.5节)以获取这些列表).此函数的精确效果如下:它首先在PARI堆栈上创建大小为长字长的内存块,并保存最终返回的块地址.如果堆栈已用完,则会显示一条消息,提示"PARI堆栈溢出",并引发错误.否则,它将设置PARI对象的类型和长度.实际上,它会填充其第一个代码字(z [0]).

请参见 https://pari. math.u-bordeaux.fr/pub/pari/manuals/2.7.6/libpari.pdf

在本文档的示例中,您可以看到创建带有两个元素的向量,调用大小为 l = 3 以获得合适的向量.实际数字矢量的第一个元素不是以索引0开头,而是以索引1开头(请参阅本PDF文档中的4.5.15节).

使用

git clone http://pari.math.u-bordeaux.fr/git/pari.git   

可以获取PARI的源代码.

在src/headers/parigen.h中,您可以在最后看到不同的类型.这是一个枚举,我们需要的类型是t_VEC.对应的整数是17.

因此,我们现在可以定义一个小函数,将一个小动物转换成GP向量,如下所示:

def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1

然后我们可以这样呼叫ellinit:

h = (0, 0, 0, 4, 6)
res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision)

要使用[0,0,0,4,6]参数对其进行测试,可以从命令行调用GP:

? ellinit([0, 0, 0, 4, 6], 1)
%1 = [0, 0, 0, 4, 6, 0, 8, 24, -16, -192, -5184, -19648, 110592/307, Vecsmall([1]), [Vecsmall([128, -1])], [0, 0, 0, 0, 0, 0, 0, 0]]

所引用的PDF文档的第441页上的示例的一个独立的小型Python程序可能看起来像这样:

from ctypes import *

pari = cdll.LoadLibrary("libpari.so")

pari.stoi.restype = POINTER(c_long)
pari.cgetg.restype = POINTER(POINTER(c_long))
pari.ellinit.restype = POINTER(POINTER(c_long))
pari.ellisdivisible.restype = c_long
pari.nfinit0.restype = POINTER(c_long)
pari.polcyclo_eval.restype = POINTER(c_long)
pari.fetch_user_var.restype = c_long
pari.pol_x.restype = POINTER(c_long)

(t_VEC, t_COL, t_MAT) = (17, 18, 19)  # incomplete
precision = c_long(38)

pari.pari_init(2 ** 19, 0)


def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1


def main():
    t = pari.pol_x(pari.fetch_user_var(bytes("t", "utf8")))
    Q = pari.pol_x(pari.fetch_user_var(bytes("Q", "utf8")))
    K = pari.nfinit0(pari.polcyclo_eval(11, t), c_long(0), precision)
    h = (0, -1, 1, 0, 0)
    res = pari.ellinit(t_vec(h), K, precision)
    P = (0, 0)
    y = pari.ellisdivisible(res, t_vec(P), pari.stoi(c_long(5)), byref(Q))

    pari.pari_printf(bytes("Q: %Ps\n", "utf8"), Q)

    print("ellisdivisible =", y)


if __name__ == '__main__':
    main()

测试

现在,我们可以调用Python程序,并将其与交互式GP程序的输出进行比较,它实际上会得到相同的结果:

Q: [Mod(-t^7 - t^6 - t^5 - t^4 + 1, t^10 + t^9 + t^8 + t^7 + t^6 + t^5 + t^4 + t^3 + t^2 + t + 1), Mod(-t^9 - 2*t^8 - 2*t^7 - 3*t^6 - 3*t^5 - 2*t^4 - 2*t^3 - t^2 - 1, t^10 + t^9 + t^8 + t^7 + t^6 + t^5 + t^4 + t^3 + t^2 + t + 1)]
ellisdivisible = 1

I would like to call PARI/GP from Python. I need to use ellisdivisible(E; P; n;{&Q}) function of PARI (see function no 3.15.35 on page 441 in this link:), so I have to pass 2 vectors or arrays (e.g, E = ellinit([0,-1,1,0,0], K);P = [0,0];), how I do that?

To call a PARI function(in C) of single argument/variable from Python (given by Thomas Baruchel), we have something like below -

import ctypes

# load the library
pari=ctypes.cdll.LoadLibrary("libpari.so")

# set the right return type of the functions
pari.stoi.restype = ctypes.POINTER(ctypes.c_long)
pari.nextprime.restype = ctypes.POINTER(ctypes.c_long)

# initialize the library 
pari.pari_init(2**19,0)

def nextprime(v):
  g = pari.nextprime(pari.stoi(ctypes.c_long(v))) # nextprime(argument) is a PARI function
  return pari.itos(g)



print( nextprime(456) )

For example I tried -

h=(0,0,0, 4,6)
pari.stoi.restype = ctypes.POINTER(ctypes.c_long*5)
pari.ellinit.restype = ctypes.POINTER(ctypes.c_long)
def ellinit(v):
  g = pari.ellinit(pari.stoi(ctypes.c_long(v)*5))
  return pari.itos(g)


print(ellinit(h))

I got below error -

   File "C:\Users\miron\Desktop\trash5\x\f.py", line 68, in <module>
    print( ellinit(h) )
  File "C:\Users\miron\Desktop\trash5\x\f.py", line 62, in ellinit
    g = pari.ellinit(pari.stoi(ctypes.c_long(v)*5))
TypeError: an integer is required (got type tuple)

How do I pass a tuple/array/vector? Thanks.

Edit: Failed attempt to get ellisdivisible(E; P; n;{&Q}) -

from ctypes import *

pari = cdll.LoadLibrary("C:\\Program Files\\Python37\\Pari64-2-11-3\\libpari.dll")

pari.stoi.restype = POINTER(c_long)
pari.cgetg.restype = POINTER(POINTER(c_long))
pari.ellinit.restype = POINTER(POINTER(c_long))
#-------------------------CHANGE 1
pari.ellisdivisible.restype = c_long
Q = pari.stoi(c_long(0))
#-------------------------
(t_VEC, t_COL, t_MAT) = (17, 18, 19)  # incomplete
precision = c_long(38)
pari.pari_init(2 ** 19, 0)
def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1
def main():
    h = (0, 0, 0, 0, 1)
    P=(0,0)
    res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision)
#---------------CHANGE 2
   # res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision).disc
    y = pari.ellisdivisible(res, t_vec(P), pari.stoi(c_long(5)), byref(Q))
    print(pari.itos(y))
#---------------
    for i in range(1, 13):
        print(pari.itos(res[i]))

if __name__ == '__main__':
    main()

The error is -

Traceback (most recent call last):
  File "C:\Users\miron\Desktop\trash5\x\ex - Copy (2).py", line 34, in <module>
    main()
  File "C:\Users\miron\Desktop\trash5\x\ex - Copy (2).py", line 28, in main
    print(pari.itos(y))
OSError: exception: access violation reading 0x0000000000000009

解决方案

Python tuples or C-arrays cannot be used directly because PARI is Using PARI/GP specific vectors where the type/length is encoded in the first element.

In section 4.4.1 Creation of PARI objects it says:

The basic function which creates a PARI object is GEN cgetg(long l, long t) l specifies the number of longwords to be allocated to the object, and t is the type of the object, in symbolic form (see Section 4.5 for the list of these). The precise effect of this function is as follows: it first creates on the PARI stack a chunk of memory of size length longwords, and saves the address of the chunk which it will in the end return. If the stack has been used up, a message to the effect that "the PARI stack overflows" is printed, and an error raised. Otherwise, it sets the type and length of the PARI object. In effect, it fills its first codeword (z[0]).

see https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.7.6/libpari.pdf

In the examples in this document you can see that to create a vector with two elements, it is called with the size l=3 to get a suitable vector. The first element of the actual number vector does not start with index 0, but with index 1 (see section 4.5.15 in this PDF document).

With

git clone http://pari.math.u-bordeaux.fr/git/pari.git   

the source code for PARI can be fetched.

There you can see the different types at the end in src/headers/parigen.h. It is an enum and the type we need is t_VEC. The corresponding integer is 17.

We can therefore now define a small function that converts a tupel into a GP vector like this:

def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1

We could then call ellinit so:

h = (0, 0, 0, 4, 6)
res = pari.ellinit(t_vec(h), pari.stoi(c_long(1)), precision)

In order to test it with your [0, 0, 0, 4, 6] parameters, you can call GP from the command line:

? ellinit([0, 0, 0, 4, 6], 1)
%1 = [0, 0, 0, 4, 6, 0, 8, 24, -16, -192, -5184, -19648, 110592/307, Vecsmall([1]), [Vecsmall([128, -1])], [0, 0, 0, 0, 0, 0, 0, 0]]

A small, self-contained Python program of the example on page 441 of the cited PDF document might look like this:

from ctypes import *

pari = cdll.LoadLibrary("libpari.so")

pari.stoi.restype = POINTER(c_long)
pari.cgetg.restype = POINTER(POINTER(c_long))
pari.ellinit.restype = POINTER(POINTER(c_long))
pari.ellisdivisible.restype = c_long
pari.nfinit0.restype = POINTER(c_long)
pari.polcyclo_eval.restype = POINTER(c_long)
pari.fetch_user_var.restype = c_long
pari.pol_x.restype = POINTER(c_long)

(t_VEC, t_COL, t_MAT) = (17, 18, 19)  # incomplete
precision = c_long(38)

pari.pari_init(2 ** 19, 0)


def t_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i in range(1, l):
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1


def main():
    t = pari.pol_x(pari.fetch_user_var(bytes("t", "utf8")))
    Q = pari.pol_x(pari.fetch_user_var(bytes("Q", "utf8")))
    K = pari.nfinit0(pari.polcyclo_eval(11, t), c_long(0), precision)
    h = (0, -1, 1, 0, 0)
    res = pari.ellinit(t_vec(h), K, precision)
    P = (0, 0)
    y = pari.ellisdivisible(res, t_vec(P), pari.stoi(c_long(5)), byref(Q))

    pari.pari_printf(bytes("Q: %Ps\n", "utf8"), Q)

    print("ellisdivisible =", y)


if __name__ == '__main__':
    main()

Test

Now we can call the Python program and, and compare with it with the output of the interactive GP program, it actually gives the same result:

Q: [Mod(-t^7 - t^6 - t^5 - t^4 + 1, t^10 + t^9 + t^8 + t^7 + t^6 + t^5 + t^4 + t^3 + t^2 + t + 1), Mod(-t^9 - 2*t^8 - 2*t^7 - 3*t^6 - 3*t^5 - 2*t^4 - 2*t^3 - t^2 - 1, t^10 + t^9 + t^8 + t^7 + t^6 + t^5 + t^4 + t^3 + t^2 + t + 1)]
ellisdivisible = 1

这篇关于如何使用Tuple/Array/Vector从Python(ctypes)调用PARI/GP?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆