Topcoder SRM554 Hard : TheBrickTower（矩阵快速幂）

John and Brus are building towers using toy bricks. They have an unlimited supply of bricks of C different colors. Each brick is a 1x1x1 cube. A tower of height X is a 2x2xX rectangular prism, built using 4X bricks.
John and Brus want their towers to look nice. A tower is nice if it has the following two properties:
1.There are at most K pairs of neighboring bricks with the same color. (Two bricks are neighboring if they share a common side.)
2.The height of the tower is between 1 and H, inclusive.
You are given the ints C and K and the long H. Return the number of nice towers, modulo 1,234,567,891.
(n≤4747,h≤474747474747474747,k≤7)

题解：
首先看数据，k特别小，可以在这上面搞事。
f[i][k]表示当前层的状态为i,相同个数有k个的个数，发现这层的颜色是什么是无所谓的，只需要知道相对的颜色状态就行了，去除了很多重复状态，最后其实只有15种状态（自行打表）。
那么对于每两个状态的转移都dfs一遍处理，最后把转移矩阵构造出来快速幂就好了。

#include<bits/stdc++.h>using namespace std;const int Mod=1234567891;const int mod[4]={1,4,16,64};int tot,lim,c,k;long long h;struct node{    int status,k,usd;    int pos[4];     node(){}    node(int status,int k,int usd):status(status),k(k),usd(usd){}};struct matrix{    long long a[130][130];    inline void init(int t){        memset(a,0,sizeof(a));        for(int i=0;i<=lim;i++)a[i][i]=t;    }    friend inline matrix operator *(const matrix &a,const matrix &b){        matrix c;c.init(0);        for(int i=0;i<=lim;i++){            for(int j=0;j<=lim;j++){                long long sum=0;                for(int k=0;k<=lim;k++){                    (sum+=1ll*a.a[i][k]*b.a[k][j]%Mod)%=Mod;                }                c.a[i][j]=sum;            }        }        return c;    }    friend inline matrix operator ^(matrix a,long long b){        matrix res;res.init(1);        for(;b;b>>=1,a=a*a)if(b&1)res=res*a;        return res;    } }e,a;vector<node>sta;inline int calcPos(int status,int pos){    if(pos!=3)status%=mod[pos+1];    status=status/mod[pos];    return status; }inline int calc(const node &t){    return (t.pos[0]==t.pos[3])+(t.pos[0]==t.pos[1])+(t.pos[1]==t.pos[2])+(t.pos[2]==t.pos[3]);}inline void dfs(int nowpos,int usd,int status){    if(nowpos==4){        node t;        t.status=status;t.usd=usd;        for(int i=0;i<4;i++)            t.pos[i]=calcPos(status,i);        t.k=calc(t);        if(t.k>k)return;        sta.push_back(t);        return;    }    for(int i=0;i<usd;++i)dfs(nowpos+1,usd,status+i*mod[nowpos]);    if(usd<c)dfs(nowpos+1,usd+1,status+usd*mod[nowpos]);}int vis[4],vt,sel[4],cnt1,cnt2,sta1,sta2,p1,p2;node t1,t2;inline int fac(int n,int t){return (t==1?n%Mod:(t==0?1:1ll*fac(n-1,t-1)*n%Mod));}inline void dfs2(int nowpos,int _diff,int _cnt){    if(nowpos==4){        for(int i=cnt1;i<=k;i++){            if(i+_cnt+cnt2>k)return;            int t=fac(c-t1.usd,_diff);            (e.a[p1+tot*i][lim]+=t)%=Mod;            (e.a[p1+tot*i][p2+tot*(i+_cnt+cnt2)]+=t)%=Mod;        }        return;    }    for(int i=0;i<nowpos;i++){        if(t2.pos[i]==t2.pos[nowpos]){            sel[nowpos]=sel[i];            dfs2(nowpos+1,_diff,_cnt+(sel[nowpos]==t1.pos[nowpos]));            return;        }    }    if(_diff+t1.usd<c){        sel[nowpos]=4;        dfs2(nowpos+1,_diff+1,_cnt);    }    for(int i=0;i<t1.usd;++i){        if(vis[i]==vt)continue;        vis[i]=vt;sel[nowpos]=i;        dfs2(nowpos+1,_diff,_cnt+(sel[nowpos]==t1.pos[nowpos]));        vis[i]=0;    }}inline void trans(int pos1,int pos2){    p1=pos1;p2=pos2;    cnt1=sta[pos1].k,cnt2=sta[pos2].k;    t1=sta[pos1],t2=sta[pos2];    if(cnt1+cnt2>k)return;    ++vt;    dfs2(0,0,0);}inline int Find(int C, int K, long long H) {    c=C,k=K,h=H;dfs(1,1,0);    tot=sta.size();lim=tot*(k+1);    e.init(0);    for(int i=0;i<tot;i++)        for(int j=0;j<tot;j++)            trans(i,j);    a.init(0);    for(int i=0;i<tot;i++){        node t3=sta[i];        a.a[0][i+t3.k*tot]=fac(c,t3.usd);        a.a[0][lim]+=fac(c,t3.usd);    }    e.a[lim][lim]=1;    a=a*(e^(H-1));     return a.a[0][lim]; }int main(){    int C,K;long long H;    scanf("%d%d%lld",&C,&K,&H);    cout<<Find(C,K,H)<<endl;}