Проблема дробного значения в Ctypes для PARI / GP - PullRequest
Новая
       68

Проблема дробного значения в Ctypes для PARI / GP

0 голосов
/ 27 марта 2020

Я написал код для сравнения решения sympy и PARI/GP, но когда я даю значение дроби D=13/12, я получаю ошибку TypeError: int expected instead of float.

Так что я изменился * От 1007 * до p1[i] = pari.stoi(c_float(numbers[i - 1])), но затем nfroots не дает вывода, обратите внимание, что я должен использовать дробь в A, B, C, D, которая может принимать $ 10 ^ 10 $ цифр после десятичной точки.

Как я могу решить эту проблему?

Код приведен ниже , чтобы загрузить файл libpari.dll, нажмите здесь - 

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

pari = cdll.LoadLibrary("libpari.dll")
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
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):
        #Changed c_long to c_float, but got no output
        p1[i] = pari.stoi(c_long(numbers[i - 1]))
    return p1


def Quartic_Comparison():
    x = Symbol('x')
    a=0;A=0;B=1;C=-7;D=13/12 #PROBLEM 1

    solution=solve(a*x**4+A*x**3+B*x**2+ C*x + D, x)
    print(solution)
    V=(A,B,C,D)
    P = pari.gtopoly(t_vec(V), c_long(-1))
    res = pari.nfroots(None, P)

    print("elements as long (only if of type t_INT): ")
    for i in range(1, pari.glength(res) + 1):        
         print(pari.itos(res[i]))
    return res               #PROBLEM 2

f=Quartic_Comparison()
print(f)

Ошибка - 

[0.158343724039430, 6.84165627596057]
Traceback (most recent call last):
  File "C:\Users\Desktop\PARI Function ellisdivisible - Copy.py", line 40, in <module>
    f=Quartic_Comparison()
  File "C:\Users\Desktop\PARI Function ellisdivisible - Copy.py", line 32, in Quartic_Comparison
    P = pari.gtopoly(t_vec(V), c_long(-1))
  File "C:\Users\Desktop\PARI Function ellisdivisible - Copy.py", line 20, in t_vec
    p1[i] = pari.stoi(c_long(numbers[i - 1]))
TypeError: int expected instead of float

1 Ответ

1 голос
/ 29 марта 2020

Система типа PARI / C очень мощная и может работать с точностью, определенной пользователем. Поэтому PARI / C необходимо использовать собственную систему типов, см., Например, Реализация типов PARI https://pari.math.u-bordeaux.fr/pub/pari/manuals/2.7.6/libpari.pdf.

Все эти внутренние типы обрабатываются как указатель на long в мире PARI / C. Не обманывайтесь этим, но тип не имеет ничего общего с долго. Возможно, его лучше всего рассматривать как индекс или дескриптор, представляющий переменную, внутреннее представление которой скрыто от вызывающей стороны.

Поэтому при переключении между PARI / C world и Python вам необходимо преобразовывать типы.

Преобразование описано, например, в разделе 4.4.6 в вышеупомянутом файле PDF.

Чтобы преобразовать double в соответствующий тип PARI (= t_REAL), необходимо вызвать функцию преобразования. dbltor.

С определением 

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

можно получить вектор PARI (t_VEC), например: 

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

Пользователь -определенная точность

Но тип Python тип double имеет ограниченную точность (ищите, например, точность с плавающей запятой в стековом потоке).

Поэтому, если вы хотите работать с определенными Точность Я бы рекомендовал использовать gdiv.

Определите его, например, так: 

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

и соответственно отрегулируйте t_vec, чтобы получить вектор этих чисел PARI: 

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

Затем вам нужно использовать realroots для вычисления корней в этом случае, см. https://pari.math.u-bordeaux.fr/dochtml/html-stable/Polynomials_and_power_series.html#polrootsreal.

Вы также можете использовать rtodbl для преобразования типа PARI t_REAL обратно в double. Но то же самое относится, так как с использованием числа с плавающей запятой вы потеряете точность. Одним из решений здесь может быть преобразование результата в строку и отображение списка со строками в выводе.

Python Программа

Автономный * Программа 1090 *, которая учитывает вышеуказанные пункты, может выглядеть следующим образом: 

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)


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


def quartic_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 in range(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

Вывод на консоль будет выглядеть следующим образом: 

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

Дополнительное примечание

Не задан вопрос, но на случай, если вы хотите избежать 13/12, вы можете преобразовать формулу из

formular

до

transformed

...