dp
fi,s f i , s 表示第 i i 层的状态为 的方案数
转移是一个矩乘的形式
用矩乘快速幂优化
#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;
const int P=1234567891;
typedef int mat[130][130];
typedef long long ll;
int cb[20010][5],num[20010],t,vis[20010],vst[10];
int a[10];
inline void judge(){
vst[1]=vst[2]=vst[3]=vst[4]=0; int ct=0;
for(int i=1;i<=4;i++){
if(!vst[a[i]]) vst[a[i]]=++ct;
a[i]=vst[a[i]];
}
int val=0;
for(int i=1;i<=4;i++) val=val*5+a[i];
if(!vis[val]){
num[t]=ct;
for(int i=1;i<=4;i++)
cb[t][i]=a[i];
t++; vis[val]=1;
}
}
int n;
mat trans;
int c,k,h,cur,matc[10],fac[20010],inv[20010],cao[20010];
void dfs(int x){
if(x>4) return judge();
for(int i=1;i<=4;i++){
a[x]=i; dfs(x+1);
}
}
inline int C(int x,int y){
if(x<y) return 0;
return 1LL*fac[x]*inv[x-y]%P;
}
vector<pair<int,int> > cc;
inline void calc(int S,int T,int nb,int tot){
for(int i=1;i<=4;i++){
if(matc[cb[T][i]]==cb[S][i]) nb--;
}
if(!nb) cur=((ll)cur+tot)%P;
}
void pfs(int x,int S,int T,int nb,int tot){
if(x>num[T]) return calc(S,T,nb,tot);
for(int i=0;i<cc.size();i++){
matc[x]=cc[i].first; cao[i]++;
pfs(x+1,S,T,nb,1LL*tot*(cc[i].second-cao[i]+1)%P);
cao[i]--;
}
}
int calc(int x,int y){
if(y==n) return 1;
if(x==n) return 0;
int sx=x/k,sy=y/k,cnt=0;
cnt+=(cb[sy][1]==cb[sy][2])+(cb[sy][3]==cb[sy][4]);
cnt+=(cb[sy][1]==cb[sy][3])+(cb[sy][2]==cb[sy][4]);
if(x%k-cnt<y%k) return 0;
cur=0; cc.clear();
for(int i=1;i<=4;i++) cc.push_back(make_pair(cb[sx][i],1));
sort(cc.begin(),cc.end()); cc.erase(unique(cc.begin(),cc.end()),cc.end());
if(c-cc.size()) cc.push_back(make_pair(5,c-cc.size()));
for(int i=0;i<cc.size();i++) cao[i]=0;
cur=0; pfs(1,x/k,y/k,x%k-cnt-y%k,1); return cur;
}
inline void Pre(){
dfs(1);
fac[0]=1; for(int i=1;i<=5000;i++) fac[i]=1LL*fac[i-1]*i%P;
inv[1]=1; for(int i=2;i<=5000;i++) inv[i]=1LL*(P-P/i)*inv[P%i]%P;
inv[0]=1; for(int i=1;i<=5000;i++) inv[i]=1LL*inv[i]*inv[i-1]%P;
}
mat res,init;
inline int count(int x){
return (cb[x][1]==cb[x][2])+(cb[x][3]==cb[x][4])+(cb[x][1]==cb[x][3])+(cb[x][2]==cb[x][4]);
}
inline void mul(mat a,mat b,mat &c){
static mat tmp;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++) tmp[i][j]=0;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++)
for(int k=0;k<=n;k++)
tmp[i][j]=(tmp[i][j]+1LL*a[i][k]*b[k][j])%P;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++) c[i][j]=tmp[i][j];
}
inline void add(mat a,mat b,mat &c){
static mat tmp;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++) tmp[i][j]=0;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++)
tmp[i][j]=((ll)a[i][j]+b[i][j])%P;
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++) c[i][j]=tmp[i][j];
}
class TheBrickTowerHardDivOne{
public:
int find(int Cc, int K, ll H){
c=Cc; k=K+1; h=H;
Pre(); n=t*(K+1);
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++)
trans[i][j]=calc(i,j);
for(int i=0;i<t;i++)
if(count(i)<=K){
//cerr<<'\t'<<i<<' '<<K-count(i)<<endl;
init[0][i*(K+1)+K-count(i)]=C(c,num[i]);
//init[1][n]=((ll)init[1][n]+C(c,num[i]))%P;
}
//for(int i=1;i<=H;i++)
// mul(init,trans,init);
for(int i=0;i<=n;i++) res[i][i]=1;
for(;H;H>>=1,mul(trans,trans,trans)) if(H&1) mul(trans,res,res);
mul(init,res,res);
return res[0][n];
}
}Main;
int main(){
cout<<Main.find(4,7,47);
for(;;);
}