Experiments homework.

homework.py

from sympy import *
from sympy.polys import subresultants_qq_zz

x_var = Symbol('x')
p_var = Symbol('p')
q_var = Symbol('q')

f_x = Matrix([
    [x_var ** 0, x_var ** 1, x_var ** 2]
])
experim_x = [-1, 0, 1]
experim_w = [p_var, 1-p_var-q_var, q_var]

m = zeros(3)
for x_i, w_i in zip(experim_x, experim_w):
    f_x_subs = f_x.subs(x_var, x_i)
    a = f_x_subs.T * f_x_subs
    m = m + w_i * a

x = 3
points = [pow(x, i) for i in range(3)]
[x_0, x_1, x_2] = points

f = Matrix(
    points
)
f_t = f.T
M = m
M_inv = M ** -1
d = f_t * M_inv * f
d_det = simplify(d.det())

G = simplify(diff(d_det, p_var))
H = simplify(diff(d_det, q_var))

G_num, _ = fraction(G)
H_num, _ = fraction(H)


H_p_coeffs = Poly(H_num, p_var).coeffs()
G_p_coeffs = Poly(G_num, p_var).coeffs()

r_matrix_p = subresultants_qq_zz.sylvester(H_num, G_num, p_var)
r_matrix_q = subresultants_qq_zz.sylvester(H_num, G_num, q_var)

r_G_H_q = r_matrix_q.det().factor()
r_G_H_p = r_matrix_p.det().factor()

p_variants = solve(r_G_H_q, p_var)
q_variants = solve(r_G_H_p, q_var)
p_c = p_variants[1]
q_c = q_variants[1]
r_c = 1 - (q_c + p_c)
D = d_det.subs([(p_var, p_c), (q_var, q_c)])

print(f'x: {x}')
print(f'M: {M}')
print(f'det f.T * invM * f = {d_det}')
print()
print(f'G: {G}')
print(f'G_num: {G_num}')
print()
print(f'H: {H}')
print(f'H_num: {H_num}')
print()
print(f'R(G, H, q) = {r_G_H_q}')
print(f'R(G, H, p) = {r_G_H_p}')
print()
print(f'p: {p_c}')
print(f'q: {q_c}')
print(f'r: {r_c}')
print(f'D: {D}')