Search code examples
pythonpari

Sending a Polynomial to PARI/GP from Python (ctypes)


I would like to call nfroots({nf}; x) function of PARI/GP from Python. (see function no 3.13.135.on page 371 in this link:), but the probllem is, I couldn't send the algebraic expression or the polynomial, that need to be send, for example, x^2-7x+12, here is a very simple example of what gp can do with a quartic polynomial:

> V = readvec("coeff.txt");
> print(V)
[1,-7,12]
> P = Pol(V);  # I get following error when I use Pol in my code:    func=self._FuncPtr((name_or_ordinal, self)) AttributeError: function 'pol' not found 
> print(P)
x^2 -7*x +12
> print(nfroots(,P))
>4, 3

From the answer of Stephan Schlecht (click here), I manage to write -

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

pari.stoi.restype = POINTER(c_long)
pari.cgetg.restype = POINTER(POINTER(c_long))

pari.nfroots.restype = POINTER(POINTER(c_long))


pari.pari_init(2 ** 19, 0)

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

def main():    
    h = "x^2-7x+12"
    res = pari.nfroots(t_vec(h))  
for i in range(1, len(res)):
         print(pari.itos(res[i]))
if __name__ == '__main__':
    main()

Note that there is specific process to create of PARI objects (see the answer of Stephan Schlecht), I changed the value for t_POL = 10 , but the code didn't work, How can I execute the above PARI/GP code from python?


Solution

  • One solution could be:

    • use gtopoly, return type is POINTER(c_long)
    • return type of nfroots is POINTER(POINTER(c_long))
    • output of result with .pari_printf

    Code

    from ctypes import *
    
    pari = cdll.LoadLibrary("libpari.so")
    
    pari.stoi.restype = POINTER(c_long)
    pari.cgetg.restype = POINTER(POINTER(c_long))
    pari.gtopoly.restype = POINTER(c_long)
    pari.nfroots.restype = POINTER(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():
        V = (1, -7, 12)
        P = pari.gtopoly(t_vec(V), c_long(-1))
        res = pari.nfroots(None, P)
        pari.pari_printf(bytes("%Ps\n", "utf8"), res)
    
    
    if __name__ == '__main__':
        main()
    

    Test

    If you run the program you get the desired output in the debug console:

    [3, 4]
    

    Conversions

    With glength one can determine the length of a vector, see https://pari.math.u-bordeaux.fr/dochtml/html/Conversions_and_similar_elementary_functions_or_commands.html#length

    With itos a long can be returned if the parameter is of type t_INT, see section 4.4.6 of https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.7.6/libpari.pdf.

    In code it would look like this:

    pari.glength.restype = c_long
    pari.itos.restype = c_long
    ... 
    print("elements as long (only if of type t_INT): ")
    for i in range(1, pari.glength(res) + 1):
        print(pari.itos(res[i]))
    

    To GENtostr gives a string representation of the argument. It could be used like so:

    pari.GENtostr.restype = c_char_p
    ...
    print("elements as generic strings: ")
    for i in range(1, pari.glength(res) + 1):
        print(pari.GENtostr(res[i]).decode("utf-8"))
    

    There are many more conversion options, see the two links above.