АОВС алгоритм Штрассена

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

#include <iostream>
#include <vector>
#include <algorithm>

using std::vector;

const int min_fast_mult_size = 64;

class matrix
{
    vector<vector<double>> m;
    matrix _get_block(int r, int c, int size)const;
    void _write_block(int r, int c, const matrix& block);
    matrix _fastmult(const matrix& arg)const;
public:
    matrix(int n) { m.assign(n, vector<double>(n)); }
    const vector<double>& operator[](int row)const { return m[row]; }
    vector<double>& operator[](int row) { return m[row]; }
    matrix operator+(const matrix& arg)const;
    matrix operator-(const matrix& arg)const;
    friend matrix mult(const matrix& arg1, const matrix& arg2); //обычное умножение
    matrix operator*(const matrix& arg)const; //быстрое умножение
    bool operator==(const matrix& arg)const { return m == arg.m; }
};

int main()
{
    matrix  A(500), B(500);
    for (size_t i = 0; i < 500; i++)
    {
        for (size_t j = 0; j < 500; j++)
        {
            A[i][j] = rand() % 10;
            B[i][j] = rand() % 10;
        }
    }
    matrix R1 = mult(A, B);
    matrix R2 = A * B;
    std::cout << (R1==R2) << "\n";
}

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

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

matrix matrix::_get_block(int r, int c, int size) const
{
    matrix block(size);
    for (size_t i = 0; i < size; i++)
        for (size_t j = 0; j < size; j++)
            block[i][j] = m[r + i][c + j];
    return block;
}

void matrix::_write_block(int r, int c, const matrix& block)
{
    int size = block.m.size();
    for (size_t i = 0; i < size; i++)
        for (size_t j = 0; j < size; j++)
            m[r + i][c + j] = block[i][j];
}

matrix matrix::_fastmult(const matrix& arg) const
{
    int n = m.size();
    int half_n = n / 2;
    if (n <= min_fast_mult_size)
        return mult(*this, arg);
    matrix a = _get_block(0, 0, half_n);
    matrix b = _get_block(0, half_n, half_n);
    matrix c = _get_block(half_n, 0, half_n);
    matrix d = _get_block(half_n, half_n, half_n);
    matrix e = arg._get_block(0, 0, half_n);
    matrix g = arg._get_block(0, half_n, half_n);
    matrix f = arg._get_block(half_n, 0, half_n);
    matrix h = arg._get_block(half_n, half_n, half_n);
    matrix P1 = a._fastmult(g - h);
    matrix P2 = (a+b)._fastmult(h);
    matrix P3 = (c+d)._fastmult(e);
    matrix P4 = d._fastmult(f - e);
    matrix P5 = (a+d)._fastmult(e + h);
    matrix P6 = (b - d)._fastmult(f + h);
    matrix P7 = (a - c)._fastmult(e + g);
    matrix res(n);
    res._write_block(0, 0, P5 + P6 + P4 - P2);
    res._write_block(0, half_n, P1 + P2);
    res._write_block(half_n, 0, P3 + P4);
    res._write_block(half_n, half_n, P5 + P1 -P7 - P3);
    return res;
}

matrix matrix::operator+(const matrix& arg) const
{
    int n = m.size();
    matrix res(n);
    for (size_t i = 0; i < n; i++)
        for (size_t j = 0; j < n; j++)
            res[i][j] = m[i][j] + arg[i][j];
    return res;
}

matrix matrix::operator-(const matrix& arg) const
{
    int n = m.size();
    matrix res(n);
    for (size_t i = 0; i < n; i++)
        for (size_t j = 0; j < n; j++)
            res[i][j] = m[i][j] - arg[i][j];
    return res;
}

matrix matrix::operator*(const matrix& arg) const
{
    int n = m.size();
    int count = 0, N=n;
    while (N >= min_fast_mult_size)
    {
        N /= 2;
        ++count;
    }
    N = (N + 1) << count;
    matrix X(N), Y(N);
    X._write_block(0, 0, *this);
    Y._write_block(0, 0, arg);
    return (X._fastmult(Y))._get_block(0,0,n);
}

matrix mult(const matrix& arg1, const matrix& arg2)
{
    int n = arg1.m.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] += arg1[i][k] * arg2[k][j];
    return res;
}