CH3401 石头游戏（矩阵快速幂加速递推）

题目链接：传送门

题目：

3401 石头游戏 0x30「数学知识」例题
描述
石头游戏在一个 n 行 m 列 (1≤n,m≤8) 的网格上进行，每个格子对应一种操作序列，操作序列至多有10种，分别用0~9这10个数字指明。
操作序列是一个长度不超过6且循环执行、每秒执行一个字符的字符串。每秒钟，所有格子同时执行各自操作序列里的下一个字符。序列中的每个字符是以下格式之一：

    数字0~9：表示拿0~9个石头到该格子。
    NWSE：表示把这个格子内所有的石头推到相邻的格子，N表示上方，W表示左方，S表示下方，E表示右方。
    D：表示拿走这个格子的所有石头。

给定每种操作序列对应的字符串，以及网格中每个格子对应的操作序列，求石头游戏进行了 t 秒之后，石头最多的格子里有多少个石头。在游戏开始时，网格是空的。
输入格式

第一行4个整数n, m, t, act。
接下来n行，每行m个字符，表示每个格子对应的操作序列。
最后act行，每行一个字符串，表示从0开始的每个操作序列。
输出格式

一个整数：游戏进行了t秒之后，所有方格中最多的格子有多少个石头。
样例输入

1 6 10 3
011112
1E
E
0

样例输出

3

样例解释

这是另一个类似于传送带的结构。左边的设备0间隔地产生石头并向东传送。设备1向右传送，直到设备2。10秒后，总共产生了5个石头，2个在传送带上，3个在最右边。

题目注：

　　答案会超int，要用longlong，看了眼数据t的范围好像是1e9，其他的无所谓了。

思路：

（以下参考李煜东《算法竞赛进阶指南》）

　　这题数据范围都没给齐，要是直接做肯定一脸懵逼。不过《算法竞赛进阶指南》上作为了矩阵快速幂的例题，那就用矩阵快速幂了。。

　　要找两个东西，状态矩阵和转移矩阵。

状态矩阵：

　　题目中的状态是一个n*m的矩阵，而一般情况下矩阵快速幂要求用一维向量作为状态矩阵，所以就把这个状态拉成长度为n*m的向量了。。。

　　原矩阵的(i, j) 为状态矩阵中的F[(i-1)*m + j]，不妨令num(i, j) = (i-1)*m + j。那么在状态矩阵中就是F[num(i, j)]了。

　　另外，令F[0] = 1，作为常数1，方便转移出常数。

转移矩阵：

　　假设ch为原矩阵对应位置i, j对应时刻k的操作：

　　① 如果ch为"N"，则Ak(num(i, j), num(i-1, j)) = 1; ch 为 "W"、"E"、"S"也是同理。被转移的位置判断一下是否合法。

　　② 如果ch为[0-9]，则A_k(0, num(i, j)) = ch-'0'（加上数值），A_k(num(i, j), num(i, j)) = 1（保留上个状态的值）

　　③ 保证F[0]不变，A_k(0, 0) = 1；

　　④ A_k的其它部分为0。

不同的时刻有不同的转移矩阵，但是由于[1-6]的最小公倍数为60（6为操作序列的最大长度），所以直接处理k = 0-59时刻的所有状态矩阵，然后分解t = q*60 + r（0 ≤ r ＜ 60），

令$A =prod_{k = 0}^{59}A_{k}$，则$F_{t} = F_{0} * A^{q} * prod_{k=0}^{r-1}A_{k}$。

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int MAX_NM = 65;

int n, m, t, act;
string opt[10];
string acts[10];
ll F[MAX_NM], A[60][MAX_NM][MAX_NM], AAA[MAX_NM][MAX_NM];

inline int num(int i, int j) {
    return (i-1)*m + j;
}

void mul(ll f[MAX_NM], ll a[MAX_NM][MAX_NM]) {
    ll c[MAX_NM];
    memset(c, 0, sizeof c);
    for (int j = 0; j < MAX_NM; j++)
        for (int k = 0; k < MAX_NM; k++)
            c[j] += f[k] * a[k][j];
    memcpy(f, c, sizeof c);
}

void mulb(ll a[MAX_NM][MAX_NM], ll b[MAX_NM][MAX_NM]) {
    ll c[MAX_NM][MAX_NM];
    memset(c, 0, sizeof c);
    for (int i = 0; i < MAX_NM; i++)
        for (int j = 0; j < MAX_NM; j++)
            for (int k = 0; k < MAX_NM; k++)
                c[i][j] += a[i][k]*b[k][j];
    memcpy(a, c, sizeof c);
}

void mulself(ll a[MAX_NM][MAX_NM]) {
    ll c[MAX_NM][MAX_NM];
    memset(c, 0, sizeof c);
    for (int i = 0; i < MAX_NM; i++)
        for (int j = 0; j < MAX_NM; j++)
            for (int k = 0; k < MAX_NM; k++)
                c[i][j] += a[i][k]*a[k][j];
    memcpy(a, c, sizeof c);
}

void init()
{
    memset(A, 0, sizeof A);
    memset(F, 0, sizeof F);
    F[0] = 1;
    for (int k = 0; k < 60; k++) {
        A[k][0][0] = 1;
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= m; j++) {
                int index = opt[i-1][j-1] - '0';
                int indey = k % acts[index].size();
                char ch = acts[index][indey];

                if (isupper(ch)) {
                    switch (ch) {
                        case 'N':
                            if (i-1 >= 1)
                                A[k][num(i, j)][num(i-1, j)] = 1; break;
                        case 'S':
                            if (i+1 <= n)
                                A[k][num(i, j)][num(i+1, j)] = 1; break;
                        case 'W':
                            if (j-1 >= 1)
                                A[k][num(i, j)][num(i, j-1)] = 1; break;
                        case 'E':
                            if (j+1 <= m)
                                A[k][num(i, j)][num(i, j+1)] = 1; break;
                        case 'D':
                            A[k][num(i, j)][num(i, j)] = 0;
                    }
                }
                if (isdigit(ch)) {
                    A[k][num(i, j)][num(i, j)] = 1;
                    A[k][0][num(i, j)] = ch - '0';
                }


            }
        }
    }
    for (int i = 0; i < MAX_NM; i++)
        AAA[i][i] = 1;
    for (int k = 0; k < 60; k++)
        mulb(AAA, A[k]);
}

int main()
{
    cin >> n >> m >> t >> act;
    for (int i = 0; i < n; i++)
        cin >> opt[i];
    for (int i = 0; i < act; i++)
        cin >> acts[i];
    init();
    int q = t/60;
    int r = t%60;
    // t = q*60 + r;
    for (; q; q >>= 1) {
        if (q&1)
            mul(F, AAA);
        mulself(AAA);
    }
    for (int i = 0; i < r; i++)
        mul(F, A[i]);
    ll ans = 0;
    for (int i = 1; i < MAX_NM; i++)
        ans = max(ans, F[i]);
    cout << ans << endl;
    return 0;
}