Возведение матрицы в степень

Question

Уже очень много обсуждали этот вопрос, но остались ещё пробелы. Написал свой код возведения матрицы в степень, и он даже работает, но в тестирующей системе падает с ошибкой выполнения на последнем тесте. Я не смог придумать такую матрицу, чтобы с ней мой код не работал, помогите найти ошибку, пожалуйста. Прикрепляю функцию перемножения двух матриц и сам код возведения. Matrix - это структура матрицы, в которой хранятся поля data, rows, cols.

Matrix Multiplication(Matrix& l, Matrix& r) {
    Matrix d;
    d.rows = r.rows;
    d.cols = l.cols;
    Allocate(d); //выделение места в памяти
    for (long long int i = 0; i < d.rows; i++)
        for (long long int j = 0; j < d.cols; j++)
        {
            d.data[i][j] = 0;
            for (int k = 0; k < r.cols; k++)
                d.data[i][j] += (l.data[i][k] * r.data[k][j]);
        }
    return d;
}

Кусок из функции main(). z - степень матрицы. array - введенная квадратная матрица.

    temp = array;
    // Если степень нулевая, возвращаем единичную матрицу
    if (z == 0) {
        for (long long int i = 0; i < temp.rows; i++) {
            for (long long int j = 0; j < temp.rows; j++) {
                if (i == j)
                    temp.data[i][j] = 1;
                else
                    temp.data[i][j] = 0;
            }
        }
    }
    else {
        while (z > 1) {
            temp = Multiplication(temp, array);
            --z;
        }
    }
    \\ Если z == 1, то возвращается исходная матрица

Полный код

#include <iostream>
#include <fstream>

//������� ���������
struct Matrix {
    long long int** data = 0;
    size_t rows = 0;
    size_t cols = 0;
};

bool Allocate(Matrix& m) {
    if (!m.data) {
        m.data = new long long int* [m.rows];
        for (size_t i = 0; i < m.rows; i++) {
            m.data[i] = new long long int[m.rows];
        }
        for (size_t y = 0; y < m.rows; y++) {
            for (size_t x = 0; x < m.rows; x++)
                m.data[x][y] = 0;
        }
        return true;
    }
    return false;
}

bool Deallocate(Matrix& m) {
    if (m.data) {
        for (size_t y = 0; y < m.rows; y++) {
            delete[] m.data[y];
        }
        delete[] m.data;
        return true;
    }
    return false;
}

void outputMatrix(std::ostream& out, const Matrix& m) {
    for (size_t y = 0; y < m.rows; y++) {
        for (size_t x = 0; x < m.cols; x++) {
            out << m.data[y][x] << ' ';
        }
        out << '\n';
    }
}

Matrix Multiplication(Matrix& lft, Matrix& rht) {
    Matrix dst;
    dst.rows = rht.rows;
    dst.cols = lft.cols;
    Allocate(dst);
    for (long long int i = 0; i < dst.rows; i++)
        for (long long int j = 0; j < dst.cols; j++)
        {
            dst.data[i][j] = 0;
            for (int k = 0; k < rht.cols; k++)
                dst.data[i][j] += (lft.data[i][k] * rht.data[k][j]);
        }
    return dst;
}


int main() {
    Matrix array;
    Matrix temp;
    long long int z, input;
    std::cin >> array.rows >> z;
    array.cols = array.rows;
    Allocate(array);
    for (long long int i = 0; i < array.rows; i++) {
        for (long long int j = 0; j < array.rows; j++) {
            std::cin >> input;
            array.data[i][j] = input;
        }
    }
    temp = array;
    while (z > 1) {
        temp = Multiplication(temp, array);
        --z;
    }
    outputMatrix(std::cout, temp);
    return 0;
}

Answer 1

Ваш код выжирает память :)

Проверял на матрице 400x400 и возводил в 100 степень

код работал 35 секунд и отъел 143МБ памяти

Написал свой, в котором

матрица была одномерной
временных матриц было минимум

код работал 6,5 секунд и отъел 5МБ памяти

потом изменил функцию, которая не умножала матрицы 100 раз, а возводила их в квадраты (на самом деле можно еще соптимизировать)

код работал 3,3 секунды и отъел 7МБ памяти

class Matrix2 {
private:
    int     m_rows;
    int     m_cols;
    int*    m_data;

public:
    Matrix2(const int width, const int height);
    Matrix2(const Matrix2& matrix);
    ~Matrix2();

    void fill(const int value);
    void identity();

    void pow(const int power);
    void pow2(const int power);
};

Matrix2::Matrix2(const int rows, const int cols) {
    m_rows = rows;
    m_cols = cols;

    const int size = m_rows * m_cols;

    m_data = new int[size];
    memset(m_data, 0, size * sizeof(int));
}

Matrix2::Matrix2(const Matrix2& matrix) {
    m_rows = matrix.m_rows;
    m_cols = matrix.m_cols;

    const int size = m_rows * m_cols;

    m_data = new int[size];
    memcpy(m_data, matrix.m_data, size * sizeof(int));
}

Matrix2::~Matrix2() {
    delete[] m_data;
}

// заполнить матрицу
void Matrix2::fill(const int value) {
    const int size = m_rows * m_cols;

    for (int index = 0; index < size; index++)
        m_data[index] = value;
}

// сделать единичную матрицу
void Matrix2::identity() {
    const int size = m_rows * m_cols;
    memset(m_data, 0, size * sizeof(int));

    for (int i = 0; i < m_rows; i++) {
        m_data[i * m_cols + i] = 1;
    }
}


// возвести матрицу в степень
void Matrix2::pow(const int power)
{
    if (power == 0)
        identity();

    if (power < 2)
        return;

    const int size = m_rows * m_cols;

    int* tmp_src = new int[size];
    int* tmp_dst = new int[size];

    memcpy(tmp_src, m_data, size * sizeof(int));

    for (int index = 1; index < power; index++) {

        for (int i = 0; i < m_rows; i++) {
            for (int j = 0; j < m_cols; j++) {

                int res = 0;
                for (int k = 0; k < m_cols; k++)
                    res += tmp_src[i * m_cols + k] * m_data[k * m_cols + j];

                tmp_dst[i * m_cols + j] = res;
            }
        }

        int* tmp = tmp_src;
        tmp_src = tmp_dst;
        tmp_dst = tmp;
    }

    memcpy(m_data, tmp_src, size * sizeof(int));

    delete[] tmp_src;
    delete[] tmp_dst;
}

void Matrix2::pow2(const int power)
{
    if (power == 0)
        identity();

    if (power < 2)
        return;

    Matrix2 matrix1(*this);
    Matrix2 matrix2(*this);

    matrix1.pow(power / 2);

    for (int i = 0; i < m_rows; i++) {
        for (int j = 0; j < m_cols; j++) {

            int res = 0;
            for (int k = 0; k < m_cols; k++)
                res += matrix1.m_data[i * m_cols + k] * matrix1.m_data[k * m_cols + j];

            m_data[i * m_cols + j] = res;
        }
    }

    if (power % 2 == 1) {

        for (int i = 0; i < m_rows; i++) {
            for (int j = 0; j < m_cols; j++) {

                int res = 0;
                for (int k = 0; k < m_cols; k++)
                    res += matrix1.m_data[i * m_cols + k] * matrix2.m_data[k * m_cols + j];

                m_data[i * m_cols + j] = res;
            }
        }
    }
}

посмотрите принцип/подход

единственное - матрицы могут быть транспонированные - код написан больше для проверки способов вычисления степеней и его надо еще причесывать :)

Answer 2

Чтоб не жрало память и работало быстро - сделаем полноценный класс, а возведение в степень сделаем быстрым:

#include <iostream>
#include <fstream>
#include <cassert>

class Matrix {
public:
    Matrix(size_t rows_ = 1, size_t cols_ = 1,
           long long int val = 0):rows_(rows_),cols_(cols_)
    {
        data = new long long int*[rows_];
        for(size_t i = 0; i < rows_; ++i)
        {
            data[i] = new long long int[cols_];
            for(size_t j = 0; j < cols_; ++j)
                data[i][j] = val;
        }
    }
    ~Matrix()
    {
        for(size_t i = 0; i < rows_; ++i)
            delete[] data[i];
        delete[] data;
    }
    Matrix(const Matrix& M):rows_(M.rows_),cols_(M.cols_)
    {
        data = new long long int*[rows_];
        for(size_t i = 0; i < rows_; ++i)
        {
            data[i] = new long long int[cols_];
            for(size_t j = 0; j < cols_; ++j)
                data[i][j] = M.data[i][j];
        }
    }

    Matrix& operator = (const Matrix& M)
    {
        Matrix tmp(M);
        swap(tmp);
        return *this;
    }

    size_t rows() const { return rows_; }
    size_t cols() const { return cols_; }

    long long * operator[](size_t R) { return data[R]; }

    friend std::ostream& operator << (std::ostream& out, const Matrix& m);
    friend Matrix operator *(const Matrix& lft, const Matrix& rht);

    Matrix pow(unsigned int z) const;

private:
    void swap(Matrix& M)
    {
        ::std::swap(data,M.data);
        ::std::swap(rows_,M.rows_);
        ::std::swap(cols_,M.cols_);
    }
    long long int** data = 0;
    size_t rows_;
    size_t cols_;
};

std::ostream& operator << (std::ostream& out, const Matrix& m)
{
    for (size_t y = 0; y < m.rows_; y++) {
        for (size_t x = 0; x < m.cols_; x++) {
            out << m.data[y][x] << ' ';
        }
        out << '\n';
    }
    return out;
}

Matrix operator *(const Matrix& lft, const Matrix& rht)
{
    assert(lft.cols_ == rht.rows_);
    Matrix dst(lft.rows_,rht.cols_,0);
    for(size_t i = 0; i < dst.rows_; i++)
        for(size_t j = 0; j < dst.cols_; j++)
        {
            for(size_t k = 0; k < rht.rows_; k++)
                dst.data[i][j] += lft.data[i][k] * rht.data[k][j];
        }
    return dst;
}

Matrix Matrix::pow(unsigned int z) const
{
    assert(rows_ == cols_);
    Matrix res(rows_,cols_), t(*this);
    for(size_t i = 0; i < rows_; ++i) res[i][i] = 1;

    for(;;)
    {
        if (z&1) res = res * t;
        z >>= 1;
        if (z) t = t*t;
        else break;
    }
    return res;
}



int main()
{
    size_t sz; unsigned int z, zz;
    std::cin >> sz >> z;
    zz = z;

    Matrix array(sz,sz);
    for(size_t i = 0; i < array.rows(); i++) {
        for(size_t j = 0; j < array.cols(); j++) {
            long long int input;
            std::cin >> input;
            array[i][j] = input;
        }
    }
    Matrix temp (array);
    temp = array;
    while (z > 1) {
        temp = temp * array;
        --z;
    }
    std::cout << temp;

    std::cout << "---------------------\n";
    temp = array.pow(zz);
    std::cout << temp;

}

БЛОГ НА HUSL

Возведение матрицы в степень

Ответы (2 шт):