[bzoj2004][Hnoi2010]Bus 公交线路——动态规划+矩阵快速幂

题目大意：

小Z所在的城市有N个公交车站，排列在一条长(N-1)km的直线上，从左到右依次编号为1到N，相邻公交车站间的距离均为1km。 作为公交车线路的规划者，小Z调查了市民的需求，决定按下述规则设计线路：

设共K辆公交车，则1到K号站作为始发站，N-K+1到N号台作为终点站。
每个车站必须被一辆且仅一辆公交车经过（始发站和终点站也算被经过）。
公交车只能从编号较小的站台驶往编号较大的站台。
一辆公交车经过的相邻两个站台间距离不得超过Pkm。
在最终设计线路之前，小Z想知道有多少种满足要求的方案。由于答案可能很大，你只需求出答案对30031取模的结果。

思路：

显然前面的k个车站和后面的k个车站需要特殊处理，然后中间的地方可以直接dp。
dp状态记录每一辆车目前离这个点最近的停靠站是哪个，转移的规则是每辆车相邻两个停靠车站之间距离\(\leq p\)。
考虑这样的状态总共有多少个，每一个点都会有一辆车停靠，于是一定会有一个0，然后每一辆车的停靠点距离这个点距离应该\(< p\)，于是最终的状态数只有\({p-1\choose k}\)，直接将所有的状态搜出来，然后处理一下状态之间转移的关系后矩阵快速幂即可。

#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
#define MREP(i,x) for(int i=beg[x],v;v=to[i],i;i=las[i])
#define debug(x) cout<<#x<<"="<<x<<endl
#define fi first
#define se second
#define mk make_pair
#define pb push_back
typedef long long ll;

using namespace std;

void File(){
	freopen("bzoj2004.in","r",stdin);
	freopen("bzoj2004.out","w",stdout);
}

template<typename T>void read(T &_){
	T __=0,mul=1; char ch=getchar();
	while(!isdigit(ch)){
		if(ch=='-')mul=-1;
		ch=getchar();
	}
	while(isdigit(ch))__=(__<<1)+(__<<3)+(ch^'0'),ch=getchar();
	_=__*mul;
}

const int mod=30031;
const int maxn=200+10;
int n,p,k,ans;
int st[maxn][11],cnt;

struct Matrix{
	int c[maxn][maxn];
	Matrix(){REP(i,1,cnt)REP(j,1,cnt)c[i][j]=0;}
};

Matrix operator * (Matrix a,Matrix b){
	Matrix ret;
	REP(i,1,cnt)REP(j,1,cnt)REP(l,1,cnt)
		ret.c[i][j]=(ret.c[i][j]+a.c[i][l]*b.c[l][j])%mod;
	return ret;
}

namespace Find_States{
	int now[11];
	void dfs(int s,int f){
		if(s>k){
			++cnt;
			REP(i,1,k)st[cnt][i]=now[i];
			return;
		}
		REP(i,f,p-1){
			now[s]=i;
			dfs(s+1,i+1);
		}
	}
}

void Match(Matrix &a){
	int t[11];
	REP(i,1,cnt){
		REP(l,1,k){
			bool lim=false;
			REP(r,1,k)if(r==l)t[r]=0;
			else{
				t[r]=st[i][r]+1;
				if(t[r]>=p){lim=true;break;}
			}
			if(lim)continue;

			sort(t+1,t+k+1);
			REP(j,1,cnt){
				bool flag=true;
				REP(r,1,k)if(t[r]!=st[j][r]){flag=false;break;}
				if(flag)a.c[i][j]=1;
			}
		}
	}
}

Matrix qpow(Matrix x,int y){
	Matrix ret;
	REP(i,1,cnt)ret.c[i][i]=1;
	while(y){
		if(y&1)ret=ret*x;
		x=x*x;
		y>>=1;
	}
	return ret;
}

void work(){
	Matrix a,t;
	t.c[1][1]=1;
	Match(a);
	a=t*qpow(a,n-k*2);

	REP(i,1,cnt){
		DREP(j,k,1)a.c[1][i]=a.c[1][i]*(min(p-st[i][j],k)-k+j)%mod;
		ans=(ans+a.c[1][i])%mod;
	}

	printf("%d\n",ans);
}

int main(){
	File();

	read(n),read(k),read(p);

	Find_States::dfs(2,1);

	work();

	return 0;
}