How to use Tuple/Array/Vector to call PARI/GP from Python (ctypes)?

Question

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

Answer 1

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