Fraction Value Problem In Ctypes To Pari/gp

I have written a code to compare the solution of sympy and PARI/GP, but when I give a fraction value D=13/12, I get error, TypeError: int expected instead of float. So I changed p

Solution 1:

The PARI/C type system is very powerful and can also work with user-defined precision. Therefore PARI/C needs to use its own types system, see e.g. Implementation of the PARI typeshttps://pari.math.u-bordeaux.fr/pub/pari/manuals/2.7.6/libpari.pdf.

All these internal types are handled as pointer to long in the PARI/C world. Don't be fooled by this, but the type has nothing to do with long. It is perhaps best thought of as an index or handle, representing a variable whose internal representation is hidden from the caller.

So whenever switching between PARI/C world and Python you need to convert types.

Conversion are described e.g. in section 4.4.6 in the above mentioned PDF file.

To convert a double to the corresponding PARI type (= t_REAL) one would therefore call the conversion function dbltor.

With the definition of

pari.dbltor.restype = POINTER(c_long)
pari.dbltor.argtypes = (c_double,)

one could get a PARI vector (t_VEC) 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.dbltor(numbers[i - 1])
    return p1

User-defined Precision

But the type Python type double has limited precision (search e.g. for floating point precision on stackoverflow).

Therefore if you want to work with defined precision I would recommend to use gdiv.

Define it e.g. like so:

V = (pari.stoi(A), pari.stoi(B), pari.stoi(C), pari.gdiv(pari.stoi(13), pari.stoi(12)))

and adjust t_vec accordingly, to get a vector of these PARI numbers:

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] = numbers[i - 1]
    return p1

You then need to use realroots to calculate the roots in this case, see https://pari.math.u-bordeaux.fr/dochtml/html-stable/Polynomials_and_power_series.html#polrootsreal.

You could likewise use rtodbl to convert a PARI type t_REAL back to a double. But the same applies, since with using a floating point number you would loose precision. One solution here could be to convert the result to a string and display the list with the strings in the output.

Python Program

A self-contained Python program that considers the above points might look like this:

from ctypes import *
from sympy.solvers import solve
from sympy import Symbol

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

pari.stoi.restype = POINTER(c_long)
pari.stoi.argtypes = (c_long,)

pari.cgetg.restype = POINTER(POINTER(c_long))
pari.cgetg.argtypes = (c_long, c_long)

pari.gtopoly.restype = POINTER(c_long)
pari.gtopoly.argtypes = (POINTER(POINTER(c_long)), c_long)

pari.dbltor.restype = POINTER(c_long)
pari.dbltor.argtypes = (c_double,)

pari.rtodbl.restype = c_double
pari.rtodbl.argtypes = (POINTER(c_long),)

pari.realroots.restype = POINTER(POINTER(c_long))
pari.realroots.argtypes = (POINTER(c_long), POINTER(POINTER(c_long)), c_long)

pari.GENtostr.restype = c_char_p
pari.GENtostr.argtypes = (POINTER(c_long),)

pari.gdiv.restype = POINTER(c_long)
pari.gdiv.argtypes = (POINTER(c_long), POINTER(c_long))

(t_VEC, t_COL, t_MAT) = (17, 18, 19)  # incomplete
precision = c_long(38)
pari.pari_init(2 ** 19, 0)


deft_vec(numbers):
    l = len(numbers) + 1
    p1 = pari.cgetg(c_long(l), c_long(t_VEC))
    for i inrange(1, l):
        p1[i] = numbers[i - 1]
    return p1


defquartic_comparison():
    x = Symbol('x')
    a = 0
    A = 0
    B = 1
    C = -7
    D = 13 / 12
    solution = solve(a * x ** 4 + A * x ** 3 + B * x ** 2 + C * x + D, x)
    print(f"sympy: {solution}")

    V = (pari.stoi(A), pari.stoi(B), pari.stoi(C), pari.gdiv(pari.stoi(13), pari.stoi(12)))
    P = pari.gtopoly(t_vec(V), -1)
    roots = pari.realroots(P, None, precision)
    res = []
    for i inrange(1, pari.glength(roots) + 1):
        res.append(pari.GENtostr(roots[i]).decode("utf-8"))  #res.append(pari.rtodbl(roots[i]))return res


f = quartic_comparison()
print(f"PARI: {f}")

Test

The output on the console would look like:

sympy: [0.158343724039430, 6.84165627596057]
PARI: ['0.15834372403942977487354358292473161327', '6.8416562759605702251264564170752683867']

Side Note

Not really asked in the question, but just in case you want to avoid 13/12 you could transform your formula from

to