C++【矩阵加速】POJ3233 Matrix Power Series

最新推荐文章于 2024-04-04 13:40:17 发布

原创最新推荐文章于 2024-04-04 13:40:17 发布 · 277 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#数论 #矩阵 #矩阵加速

My OI 专栏收录该内容

32 篇文章

订阅专栏

该博客介绍了如何利用矩阵加速技术解决POJ3233问题，通过矩阵乘法简化求解Sk=∑i=1kximodm的过程。难点包括矩阵套矩阵的操作及单位矩阵和零矩阵的应用。

文章目录

题目描述
前置芝士
问题的简化
本题做法
代码（你看不到我）
附：头文件集锦

题目描述

POJ-3233
看到这令人目瞪口呆的数据范围， $O(\log n)$ 没逃了，于是，最好想到的肯定是矩阵加速。

前置芝士

矩阵，矩阵乘法，矩阵加速

~~到这里找就好了嘛~~

问题的简化

首先抛掉矩阵，先想一想这样一道题：
给定 $m, k, x$ ，求 $S_k=\sum_{i=1}^{k}x^i\bmod m$ 。

这道题可以轻松想出两种对于矩阵加速状态的定义： $\begin{bmatrix}S_i&x^i\end{bmatrix}$ 或 $\begin{bmatrix}S_i&1\end{bmatrix}$ 。
对于两种状态~~，都很容易构造出~~的加速矩阵，如下：
$\begin{bmatrix}1&0\\x&x\end{bmatrix}\qquad and\qquad\begin{bmatrix}x&0\\x&1\end{bmatrix}$
于是这道题就AC_？了……

本题做法

就和上面的题一样，就是把原先的 $x$ 变成了矩阵 $A$ 。

难点1

原先的 $x$ 是要放在一个矩阵中的，那把矩阵 $A$ 还要套在矩阵中吗？怎么套？
这是可以实现的，使用矩阵套矩阵的方法，模板代码如下：

int N, M;
struct Matrix{
    int n, m;
    long long A[31][31];
    inline Matrix(int N = 0, int M = 0) {
        n = N; m = M; memset(A, 0, sizeof A);
    }
    inline Matrix operator * (const Matrix &rhs) const {
        Matrix ans = Matrix(n, rhs.m);
        for(reg int k = 1; k <= m; k++)
            for(reg int i = 1; i <= ans.n; i++)
                for(reg int j = 1; j <= ans.m; j++)//据说kij比ijk快（Floyd？）
                    ans.A[i][j] = ans.A[i][j] + A[i][k] * rhs.A[k][j];
        return ans;
    }
    inline Matrix operator + (const Matrix & rhs) const {
        Matrix ans = Matrix(n, m);
        for(reg int i = 1; i <= n; i++)
            for(reg int j = 1; j <= m; j++)
                ans.A[i][j] = (A[i][j] + rhs.A[i][j]) % MOD;
        return ans;
    }
};
struct Mmatrix{
    int n, m;
    Matrix A[3][3];
    inline Mmatrix(int N = 0, int M = 0) {
        n = N; m = M; memset(A, 0, sizeof A);
    }
    inline Mmatrix operator * (const Mmatrix& rhs) const {
        Mmatrix ans = Mmatrix(n, rhs.m);
        for(int k = 1; k <= m; k++)
            for(int i = 1; i <= ans.n; i++)
                for(int j = 1; j <= ans.m; j++) {
                    if(!ans.A[i][j].n) ans.A[i][j] = Matrix(N, M);
                    //注意这一步很关键，因为这个时候ans.A[i][j]还是一个0行0列的矩阵
                    //需要把它赋上行和列的值，才能进行加法运算，否则答案就会为0
                    ans.A[i][j] = ans.A[i][j] + A[i][k] * rhs.A[k][j];
                }
        return ans;
    }
};

难点2

原先定义的矩阵中是有一个 $1$ 和 $0$ 的，但是到了矩阵中 $1$ 和 $0$ 要变成什么呢？

因为矩阵加速用的是乘法运算，所以令表示 $1$ 的矩阵为 $E$ ，则 $E$ 要做到的就是 $E\times A$ （ $A$ 为任意矩阵） $= A$ 。
所以可以~~轻松的~~得到 $E$ 是一个对角线的矩阵，又称单位矩阵。
至于 $0$ 的代替物，就是全为 $0$ 的矩阵辣！

代码（你看不到我）

方法1：~~自动省略input和output~~

int main() {
    ios::sync_with_stdio(false);
    cin.tie(); cout.tie();
    cin >> N >> k >> MOD;
    A.input(N, N);
    E = Matrix(N, N);
    for(reg int i = 1; i <= N; i++) E.A[i][i] = 1;
    O = Matrix(N, N);
 
    f = Mmatrix(1, 2);//S_1, A^1
    f.A[1][1] = A; f.A[1][2] = A;
    h = Mmatrix(2, 2);
    h.A[1][1] = A; h.A[1][2] = O;
    h.A[2][1] = E; h.A[2][2] = E;
 
    piow(k - 1);
    f.A[1][1].output();
}

方法2：

int main() {
    ios::sync_with_stdio(false);
    cin.tie(); cout.tie();
    cin >> N >> k >> MOD;
    A.input(N, N);
    E = Matrix(N, N);
    for(reg int i = 1; i <= N; i++) E.A[i][i] = 1;
    O = Matrix(N, N);

    f = Mmatrix(1, 2);
    f.A[1][1] = A; f.A[1][2] = E;
    h = Mmatrix(2, 2);
    h.A[1][1] = A; h.A[1][2] = O;
    h.A[2][1] = A; h.A[2][2] = E;

    piow(k - 1);
    f.A[1][1].output();
}

附：头文件集锦

不打这么多头文件我看不起你！！！

#include<map>
#include<set>
#include<list>
#include<queue>
#include<deque>
#include<stack>
#include<ctime>
#include<vector>
#include<bitset>
#include<cctype>
#include<string>
#include<sstream>
#include<fstream>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<iomanip>
#include<iostream>
#include<algorithm>
#include<tr1/unordered_map>
#include<tr1/unordered_set>
using namespace std;
using namespace tr1;