矩阵模板

已通过洛谷模板:(【矩阵快速幂】【行列式计算】)

template<class T>
struct Matrix {
    i64 N;
    vector<vector<T>> A;
    Matrix() : N(0) {}
    Matrix(int n) : N(n), A(n, vector<T>(n, 0)) {}
    Matrix(const vector<vector<T>>& mat) {
        assert(!mat.empty() && !mat[0].empty()); // 确保非空
        N = mat.size();
        A = mat;
    }

    const vector<T>& operator[](int i) const { return A[i]; }

    Matrix operator*(const Matrix &b) const {
        assert(N == b.N);
        Matrix res(N);

        for (int i = 0; i < N; i++) {
            for (int k = 0; k < N; k++) {
                T tmp = A[i][k];
                if (tmp == 0) {
                    continue;
                }
                for (int j = 0; j < N; j++) {
                    res.A[i][j] += tmp * b.A[k][j];
                }
            }
        }
        return res;
    }

    //快速幂实现
    Matrix power(Matrix res, i64 k) const {
        Matrix base = *this;
        while (k) {
            if (k & 1) {
                res = res * base;
            }
            base = base * base;
            k >>= 1;
        }
        return res;
    }

    //行列式计算
    T det() {
        vector<vector<T>> mat = A; // 复制矩阵
        T res = 1;
        int sign = 1; // 记录行交换的符号

        for (int i = 0; i < N; i++) {
            // 寻找主元(非零元素)
            int pivot = -1;
            for (int j = i; j < N; j++) {
                if (mat[j][i] != 0) {
                    pivot = j;
                    break;
                }
            }

            // 如果该列全零,行列式为0
            if (pivot == -1) {
                return 0;
            }

            // 交换行
            if (pivot != i) {
                swap(mat[i], mat[pivot]);
                sign = -sign;
            }

            // 使用当前行消去下方行
            for (int j = i + 1; j < N; j++) {
                // 辗转相除消元
                while (mat[j][i] != 0) {
                    i64 a = mat[i][i].val(), b = mat[j][i].val();
                    i64 ratio = a / b;
                    for (int k = i; k < N; k++) {
                        mat[i][k] -= ratio * mat[j][k];
                    }
                    swap(mat[i], mat[j]);
                    sign = -sign;
                }
            }

            // 累乘主对角线元素
            res *= mat[i][i];
        }

        res *= sign;
        return res;
    }
};
posted @ 2023-12-05 15:20  Ke_scholar  阅读(38)  评论(2)    收藏  举报