矩阵加速递推 入门指南

首先是一个例题，输入N和M，问用1X2的方块填满4XN的格子共有多少种方案数，答案对M取模。

分析问题可以找出所有状态，用一个1X4的位表示一行的方块情况，0表示暂时为空，1表示恰好填入方块。

对每一行进行递推，当上一行为0时，则下一行必须放入一个竖着的方块，置为1；

当上一行为1时，下一行可以为竖格子留空置0也可以横放置1。

虽然4个二进制位共有16中可能，但是合法的方案数实际只有6中，0000、0011、0110、1001、1100、1111。

可以将其编号压缩为0~5，不过16种方案也不是很多，不压缩也可以。

接下来推导公式，若上一行为0000，则下一行必须全部放入竖方块置为1111；

若上一行为1111，则下一行可以为0000、0011、1100、0110、1111等等。

以此类推出所有的状态，这一步可以通过两行枚举与一次判断做到，由于本题列数只有4，因此手推也是可以的。

令F[i][j]表示状态i能否到达状态j，能置为1，不能置为0。

可以发现F是一个NXN的矩阵（N是指状态数），令D(S)表示当前行状态为S的方案数，D'(S)为下一行状态为S的方案数。

那么D就是一个1XN的矩阵，由矩阵乘法的性质显然有 F X D = D'

如果将D初始化为第0行的状态的话，那么做几次乘法运算就能得到第几行的状态。

利用矩阵快速幂来加速乘法运算。

得到最后一行的状态D，那么最后一行全部填满的状态1111的方案数就是所求的答案。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;
const int maxn=1<<4;
int f[maxn][maxn];
int MOD;
void init(){
    memset(f,0,sizeof(f));
    f[0][15]=1;

    f[3][12]=1;
    f[3][15]=1;

    f[6][9]=1;
    f[6][15]=1;

    f[9][6]=1;

    f[12][3]=1;
    f[12][15]=1;

    f[15][0]=1;
    f[15][3]=1;
    f[15][6]=1;
    f[15][12]=1;
    f[15][15]=1;
}
int n,m;

struct Matrix{
    int n,m;
    int a[maxn][maxn];
    Matrix(){}
    Matrix(int _n,int _m){
        clear();
        n=_n;
        m=_m;
    }
    void clear(){
        n=m=0;
        memset(a,0,sizeof(a));
    }
    Matrix getI(int n) const{
        Matrix res(n,n);
        for (int i=0;i<n;i++) res.a[i][i]=1;
        return res;
    }
    Matrix& operator=(const Matrix& rhs){
        clear();
        n=rhs.n;
        m=rhs.m;
        for (int i=0;i<n;i++){
            for (int j=0;j<m;j++){
                a[i][j]=rhs.a[i][j];
            }
        }
        return *this;
    }
    Matrix operator*(const Matrix &b) const{
        Matrix tmp;
        tmp.clear();
        tmp.n=n;
        tmp.m=b.m;
        for (int i=0;i<n;i++){
            for (int k=0;k<m;k++){
                for (int j=0;j<b.m;j++){
                    tmp.a[i][j]+=((long long)a[i][k]*b.a[k][j])%MOD;
                    tmp.a[i][j]%=MOD;
                }
            }
        }
        return tmp;
    }
    Matrix operator^(int n) const{
        Matrix ans=getI(m);
        Matrix A=(*this);
        while (n){
            if (n&1) ans=ans*A;
            A=A*A;
            n>>=1;
        }
        return ans;
    }
};

int main()
{
    init();
    while (~scanf("%d%d",&n,&m)){
        if (n==0||m==0) break;
        MOD=m;
        Matrix A(maxn,maxn);
        Matrix F(maxn,1);
        for (int i=0;i<maxn;i++){
            for (int j=0;j<maxn;j++){
                A.a[i][j]=f[i][j];
            }
        }
        F.a[(1<<4)-1][0]=1;
        Matrix ans=(A^n)*F;
        printf("%d\n",ans.a[(1<<4)-1][0]);
    }
    return 0;
}