[HDOJ2604]Queuing（递推，矩阵快速幂）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=2604

递推式是百度的，主要是练习一下如何使用矩阵快速幂优化。

递推式：f(n)=f(n-1)+f(n-3)+f(n-4)，其中f(0)=2, f(1)=4, f(2)=6, f(3)=9。

当n>4时候，需要通过这个关系来递推。

构造矩阵这种东西我以前一直认为是很玄学的，但是如果深入研究的话不难发现其实也有规律可循。这是一个齐次递推式，很好构造。

我们希望通过如下矩阵(1)得到矩阵(2)

| f(n) |      |f(n+1)|
|f(n-1)|      | f(n) |
|f(n-2)|      |f(n-1)|
|f(n-3)|      |f(n-2)|
  (1)            (2) 

那么实际上，我们是通过一个转换矩阵把f(n),f(n-1),f(n-2)移动了一下位置，并且推出了f(n+1)。再考虑到矩阵乘法的特性，我们就能得到一个递推矩阵：

|1 0 1 1|
|1 0 0 0|
|0 1 0 0|
|0 0 1 0|

我把第一行列为递推方程的系数，这样递推矩阵的第一列和目标矩阵的整列相乘即可得到递推的下一个结果。那剩下的怎么办呢？向下看会发现，按照行x列的特性，分别又列到了对应的位置。再根据矩阵乘法的性质就会很容易证明满足递推方程。

代码如下：

#include <algorithm>
#include <iostream>
#include <iomanip>
#include <cstring>
#include <climits>
#include <complex>
#include <fstream>
#include <cassert>
#include <cstdio>
#include <bitset>
#include <vector>
#include <deque>
#include <queue>
#include <stack>
#include <ctime>
#include <set>
#include <map>
#include <cmath>

using namespace std;

const int maxn = 10;
typedef struct Matrix {
    int m[maxn][maxn];
    int r;
    int c;
    Matrix(){
        r = c = 0;
        memset(m, 0, sizeof(m));
    } 
} Matrix;

Matrix mul(Matrix m1, Matrix m2, int mod) {
    Matrix ans = Matrix();
    ans.r = m1.r;
    ans.c = m2.c;
    for(int i = 1; i <= m1.r; i++) {
        for(int j = 1; j <= m2.r; j++) {
               for(int k = 1; k <= m2.c; k++) {
                if(m2.m[j][k] == 0) continue;
                ans.m[i][k] = ((ans.m[i][k] + m1.m[i][j] * m2.m[j][k] % mod) % mod) % mod;
            }
        }
    }
    return ans;
}

Matrix quickmul(Matrix m, int n, int mod) {
    Matrix ans = Matrix();
    for(int i = 1; i <= m.r; i++) {
        ans.m[i][i]  = 1;
    }
    ans.r = m.r;
    ans.c = m.c;
    while(n) {
        if(n & 1) {
            ans = mul(m, ans, mod);
        }
        m = mul(m, m, mod);
        n >>= 1;
    }
    return ans;
}

int n, m;

int main() {
    // freopen("in", "r", stdin);
    while(~scanf("%d %d", &n, &m)) {
        Matrix p;
        p.r = 4, p.c = 1;
        p.m[1][1] = 9;
        p.m[2][1] = 6;
        p.m[3][1] = 4;
        p.m[4][1] = 2;
        if(n <= 4) {
            printf("%d
", p.m[4-n+1][1] % m);
            continue;
        }
        Matrix s;
        s.r = s.c = 4;
        s.m[1][1] = 1; s.m[1][2] = 0, s.m[1][3] = 1, s.m[1][4] = 1;
        s.m[2][1] = 1; s.m[2][2] = 0, s.m[2][3] = 0, s.m[2][4] = 0;
        s.m[3][1] = 0; s.m[3][2] = 1, s.m[3][3] = 0, s.m[3][4] = 0;
        s.m[4][1] = 0; s.m[4][2] = 0, s.m[4][3] = 1, s.m[4][4] = 0;
        s = quickmul(s, n-4, m);
        int ans = 0;
        for(int i = 1; i <= p.r; i++) {
            ans = (ans + (p.m[i][1] * s.m[1][i]) % m) % m;
        }
        printf("%d
", ans % m);
    }
    return 0;
}