Как получить собственные значения матрицы 4x4 используя QR алгоритм?

У меня появилась задача нахождения корней уравнения 4 степени. Я прочитал, что для этого в программах используют QR алгоритм нахождения собственных значений матриц-компаньонов, в частности такую систему использует numpy в функции numpy.roots(). Однако я пишу код для arduino, а там нет таких математических библиотек, способных сразу выдать все корни, поэтому я решил самостоятельно изучить эту тему и написать код на основе полученной информации. Решил написать код пока что на python так как он мне более знаком и привычен для работы, после чего я просто перезапишу его на arduino. После QR алгоритма получается новая матрица, где по идее на главной диагонали располагаются собственные значения, а то есть и сами корни уравнения, но при наличии комплексных значений алгоритм сходится не к верхней треугольной матрице, а к блочной верхней треугольной, где в блоках 2х2 располагаются комплексные значения в виде собственных значений этих блоков. Вопрос заключается в том как мне "вытянуть" эти собственные значения из получившейся матрицы и как понять когда нужно остановить итерации QR алгоритма?

import math
import numpy as np

def quadratic(a, b, c):
    D = b**2-4*a*c
    if D == 0:
        x = -b/(2*a)
        return x, x
    elif D > 0:
        x1 = (-b-math.sqrt(D))/(2*a)
        x2 = (-b+math.sqrt(D))/(2*a)
        return x1, x2
    else:
        x1 = (-b-np.sqrt(D+0j))/(2*a)
        x2 = (-b+np.sqrt(D+0j))/(2*a)
        return x1, x2

def min_len(a):
    return min([len(x) for x in a])

def my_zip(*iterables):
    l = min_len(iterables)
    res = [[] for _ in range(l)]
    for it in iterables:
        for i in range(l):
            res[i].append(it[i])
    return res

def mulyiply(*args):
    if len(args) == 1 and type(args[0]) == list or type(args[0]) == tuple:
        args = args[0]
    res = 1
    for arg in args:
        res *= arg
    return res

def sub(*args):
    if len(args) == 1 and type(args[0]) == list or type(args[0]) == tuple:
        args = args[0]
    res = args[0]
    for arg in args[1:]:
        res -= arg
    return res

def vector_projection(a, b):
    dot_ab = sum(x[0] * y[0] for x, y in my_zip(a, b))
    dot_bb = sum(x[0] * x[0] for x in b)
    scalar = dot_ab/dot_bb
    
    return [[scalar * x[0]] for x in b]

def dot(*vectors):
    return sum(mulyiply(transpose(el)) for el in transpose(list(vectors)))

def vector_sum(*vectors):
    return [[sum(to_geometry_vector(el))] for el in transpose(list(vectors))]

def vector_sub(*vectors):
    return [[sub(to_geometry_vector(el))] for el in transpose(list(vectors))]

def vector_length(vector):
    return math.sqrt(sum(x[0]**2 for x in vector))

def normalise_vector(vector):
    return [[x[0]/vector_length(vector)] for x in vector]

def col(a, indx):
    return [row[indx] for row in a]

def transpose(matrix):
    return [col(matrix, i) for i in range(len(matrix[0]))]

def to_geometry_vector(vector):
    if len(vector[0]) == 1:
        return [el[0] for el in vector]
    elif len(vector[0]) > 1:
        return vector[0]

def to_alg_vector(vector, type):
    if type == 1:
        return my_zip(vector)
    elif type == 2:
        return [vector]

def matrix_multiple(a, b):
    res = []
    b = transpose(b)
    if type(b[0]) != list:
        b = [b]
    for i in range(len(col(a, 0))):
        res.append([])
        for j in range(len(b)):
            res[i].append(sum(x*y for x, y in my_zip(a[i], b[j])))
    return res

def Gram_Schmidt_ortogonalization(F):
    G1 = []
    G = []
    for i in range(len(F[0])):
        f = my_zip(col(F, i))
        g = f
        for j in range(i):
            ng = G1[i-(j+1)]
            scalar = vector_projection(f, ng)
            g = vector_sub(g, scalar)
        G1.append(g)
        G.append(to_geometry_vector(normalise_vector(g)))
    G = transpose(G)
    return G

def QR_decompozition(A):
    Q = Gram_Schmidt_ortogonalization(A)
    R = matrix_multiple(transpose(Q), A)
    for i in range(1, len(R)):
        for j in range(i):
            R[i][j] = 0
    return Q, R

def QR_algorithm(A, max_iter=50000):
    for _ in range(max_iter):
        Q, R = QR_decompozition(A)
        A = matrix_multiple(R, Q)
    return A

Вот пример работы алгоритма с уравнением x^4-5x^3+15x^2-63x+108=0, где корнями являются числа

• 3
• -0.5+3.4278273i
• -0.5-3.4278273i

Код:

a = [[5, 15, -63, 108],
     [-1, 0, 0, 0],
     [0, 1, 0, 0],
     [0, 0, 1, 0]]

print(QR_algorithm(a, 1000))

Вывод:

>>> 
[[ 2.561813756616427,     11.050856678961864,    -31.99763823763548,  113.71173687663472],
[-1.9115896707285687,    -3.5618137566160413,    11.090746033873321, -39.439140896586785],
[                0.0, 1.0007671476289143e-59,    3.0029946167608914,  -9.503500836361985], 
[                0.0,                    0.0, 9.436236899288885e-07,   2.997005383213664]]

Матрица a уже приведена к форме Хессенберга