АОВС Вычисление корней многочлена

// ConsoleApplication3.cpp : Этот файл содержит функцию "main". Здесь начинается и заканчивается выполнение программы.
//https://pastebin.com/W4B9bPhg

#include <iostream>
#include <vector>
#include <complex>
#include <algorithm>
#include <random>
#include <numeric>
#include <iomanip>

using cdouble = std::complex<long double>;
using std::vector;

class matrix
{
    vector<vector<cdouble>> m;
public:
    matrix(size_t n) { m.assign(n, vector<cdouble>(n)); }
    size_t size() const { return m.size(); }
    vector<cdouble>& operator[](size_t row) { return m[row]; }
    const vector<cdouble>& operator[](size_t row)const { return m[row]; }
    matrix operator*(const matrix& other)const
    {
        auto n = size();
        matrix res(n);
        for (size_t i = 0; i < n; i++)
            for (size_t j = 0; j < n; j++)
                for (size_t k = 0; k < n; k++)
                    res[i][j] += m[i][k] * other[k][j];
        return res;
    }
    matrix operator^(size_t n)const
    {
        matrix res(size()), p = *this;
        for (size_t i = 0; i < size(); i++)
            res[i][i] = 1;
        while (n > 0)
        {
            if (n & 1) res = res * p;
            n /= 2;
            if (n > 0) p = p * p;
        }
        return res;
    }
};

class polynom
{
    vector<cdouble> c;
public:
    polynom(int deg = 0) { c.assign(deg + 1, 0); }
    polynom(const vector<cdouble>& coefs) { c = coefs; }
    polynom operator+(const polynom& arg)const;
    polynom operator-(const polynom& arg)const;
    polynom operator*(const polynom& arg)const;
    cdouble operator[](size_t n)const;
    cdouble& operator[](size_t n);
    cdouble operator()(cdouble arg)const; //вычисление значения многочлена в т. arg
    cdouble derivate(cdouble arg) const; //вычисление производной многочлена в т.arg
    cdouble get_root()const;
    vector<cdouble> roots()const;
    polynom shift(cdouble x_0)const; // P(x) --> P(x+x_0)
    polynom divide(cdouble r)const; // P(x) --> P(x)/(x-r)
    size_t degree()const { return c.size() - 1; }
};



int main()
{
    polynom P({ -6,2,-6,3,-3,1 });
    polynom Q({ -6,2,-9,3,-3,1 }); //(x^2+1)(x^2+2)(x-3)
    polynom R({ 7,6,5,4,3,2,1 });
    std::cout << std::setprecision(20);
    std::cout << P.get_root() << std::endl;
    std::cout << Q.get_root() << std::endl;
    std::cout << R.get_root() << std::endl;
    std::cout << sizeof(long double) << std::endl;//"Hello World!\n";
}

// Запуск программы: CTRL+F5 или меню "Отладка" > "Запуск без отладки"
// Отладка программы: F5 или меню "Отладка" > "Запустить отладку"

// Советы по началу работы
//   1. В окне обозревателя решений можно добавлять файлы и управлять ими.
//   2. В окне Team Explorer можно подключиться к системе управления версиями.
//   3. В окне "Выходные данные" можно просматривать выходные данные сборки и другие сообщения.
//   4. В окне "Список ошибок" можно просматривать ошибки.
//   5. Последовательно выберите пункты меню "Проект" > "Добавить новый элемент", чтобы создать файлы кода, или "Проект" > "Добавить существующий элемент", чтобы добавить в проект существующие файлы кода.
//   6. Чтобы снова открыть этот проект позже, выберите пункты меню "Файл" > "Открыть" > "Проект" и выберите SLN-файл.

polynom polynom::operator+(const polynom& arg) const
{
    auto deg = std::max(degree(), arg.degree());
    polynom res(deg);
    for (size_t i = 0; i <= deg; i++)
    {
        res[i] = (*this)[i] + arg[i];
    }
    return res;
}

polynom polynom::operator-(const polynom& arg) const
{
    auto deg = std::max(degree(), arg.degree());
    polynom res(deg);
    for (size_t i = 0; i <= deg; i++)
    {
        res[i] = (*this)[i] - arg[i];
    }
    return res;
}

polynom polynom::operator*(const polynom& arg) const
{
    polynom res(degree() + arg.degree());
    for (size_t i = 0; i <= degree(); i++)
    {
        for (size_t j = 0; j <= arg.degree(); j++)
        {
            res[i + j] += (*this)[i] * arg[j];
        }
    }
    return res;
}

cdouble polynom::operator[](size_t n) const
{
    if (n > degree()) return 0;
    return c[n];
}

cdouble& polynom::operator[](size_t n)
{
    return c[n];
}

cdouble polynom::operator()(cdouble arg) const
{
    int n = degree();
    auto res = c[n];
    for (int i = n - 1; i >= 0; i--)
        res = res * arg + c[i];
    return res;
}

cdouble polynom::derivate(cdouble arg) const
{
    int n = degree();
    cdouble res = c[n]*cdouble(n);
    for (int i = n - 1; i > 0; i--)
        res = res * arg + c[i]*cdouble(i);
    return res;
}

cdouble polynom::get_root() const
{
    size_t n = degree();
    std::random_device rd;
    std::mt19937_64 mtrand(rd());
    std::uniform_real_distribution<> dist(0.1, 0.5);
    cdouble a{ dist(mtrand),dist(mtrand) };
    auto P = shift(a);
    matrix A(n);
    for (size_t i = 0; i < n; i++)
    {
        A[0][n - i - 1] = -P[i] / P[n];
        if (i + 1 < n) A[i + 1][i] = 1;
    }
    for (size_t k = 0; k < 10; k++)
    {
        A = A * A;
        cdouble maxA = 0;
        for (size_t i = 0; i < n; i++)
            for (size_t j = 0; j < n; j++)
                if (abs(A[i][j]) > abs(maxA))
                    maxA = A[i][j];
        for (size_t i = 0; i < n; i++)
            for (size_t j = 0; j < n; j++)
                A[i][j] /= maxA;
    }
    std::uniform_real_distribution<> dst(0.5, 1);
    std::vector<cdouble> init_vec(n);
    for (size_t i = 0; i < n; i++)
        init_vec[i] = dst(mtrand);
    cdouble num, denom;
    num = std::inner_product(begin(A[0]), end(A[0]), begin(init_vec), cdouble(0.0));
    denom = std::inner_product(begin(A[1]), end(A[1]), begin(init_vec), cdouble(0.0));
    cdouble r= num / denom + a;
    for (int i = 0; i < 5; ++i)
        r = r - (*this)(r) / derivate(r);
    return r;
}

vector<cdouble> polynom::roots() const
{
    return vector<cdouble>();
}

polynom polynom::shift(cdouble x_0) const
{
    polynom res(0), m({ x_0,1 });
    res[0] = c[degree()];
    for (int i = (int)degree() - 1; i >= 0; i--)
    {
        res = res * m;
        res[0] += c[i];
    }
    return res;
}

polynom polynom::divide(cdouble r) const
{
    polynom res(degree() - 1);
    res.c.back() = c.back();
    for (int i = (int)degree() - 2; i >= 0; --i)
        res[i] = res[i + 1] * r + c[i];
    return res;
}