题目类型：状压(DP) -> 矩阵乘法
绝妙然而思维难度极其大的一道好题！
传送门：>Here<
题意：有一个环形花圃，可以种两种花：0或1. 要求任意相邻的(M)个花中1的个数不超过(k)个。总共有(N)个花。问方案数
解题思路
非常巧妙的一道题。
先看如何拿到(80pts)
(N leq 10^5)，也就是说可以(O(n))带若干常数。我们发现影响当前状态的决策的仅仅就是离它最近的那(M)个花圃。由此可以进行状压(DP)，(dp[i][s])表示目前决策第(i)个花圃，并且左侧(M)个花圃的状态为(s)时的方案数。很明显可以通过(i-1)时的状态来进行转移。由于仅仅只是向右移动了一格，所以原先的右边(M-1)个不动，新加进来的那个最右侧的可以是1或者0。当然在写方程的时候是要倒过来的，于是我们可以得到方程$$dp[i][s] = dp[i-1][(s/2)+2^{M-1}] + dp[i-1][s/2]$$我们可以先(dfs)预处理出所有的可能状态。
那么题目说花圃是个环形，怎么办的？其实好办。我们令(dp[M][s]=1)，然后一路转移到(dp[N+M])，这时取(dp[N+M][s])作为(s)为初始状态（前(M)个花圃）的答案。因为前(M)个花圃等同于(N+1..M)的花圃。他们的状态吻合了（都是(s)）就对了。所以我们需要所有可行的枚举(s)作为初始状态。答案累积
仅仅就是递推？那是否可以，矩阵乘法？？
矩阵乘法优化递推，然而这道题还略微有些复杂。
首先我们可以改写一下方程，使得它更加具备(Floyd)的外貌。
不如浪费一层循环，去扫一个状态(k)。使得如果(k)可以转移到(j)，那么$$dp[i][j]=sumlimits_{}dp[i-1][k]$$或者进一步，我们连(if)语句也省略掉，预处理一个布尔数组，其中(b[k][j])表示(k)能否转移到(j)。那么$$dp[i][j]=sumlimits_{}dp[i-1][k]*b[k][j]$$这样一来，这个式子就是标准的矩阵乘法了。我们可以忽略(i)的存在，它的结果就等同于初始的(dp)数组乘以这个布尔数组的(N)次方
于是我们就可以利用矩阵乘法将(b)数组做一个矩阵快速幂。
答案究竟是什么？
答案究竟是什么？我们究竟用什么来乘这个(N)次方的结果？
我们还是参考刚才的状压做法。枚举一个(s)，此时只有(dp[M][s])为1，其他都为0. 那么我们可以想象我们其实有(32)个初始的项，分别是(dp[M][0],dp[M][1],..,dp[M][s],..,dp[M][32])。其中只有(dp[M][s])为1. 每乘一次矩阵就刷新一遍，成为(dp[M+1][0],dp[M+1][1],..,dp[M+1][s],..,dp[M+1][32]).直到刷新(N)遍以后变成(dp[M+N][0],dp[M+N][1],..,dp[M+N][s],..,dp[M+N][32])。然而对于每一个(s)，我们只需要取(dp[N+M][s])作为答案。
考虑我们的(dp[N+M][s])是怎么得来的？是(dp)矩阵依次乘以布尔矩阵的第(s)列的和。而我们刚才说了，只有(dp[N+M][s])为1，其他都是0。因此结果就等同于是(b[s][s])这一项。
我们发现根本不用枚举(s)！也就是说，只需要做一遍矩阵乘法——因为对于状态(s)，对应的答案一定是(b[s][s])。那么最终的答案也就是所有的(b[s][s])？即矩阵的对角线之和！当然，并不是整一条对角线，要确保(s)是合法的。
反思
思维难度及其深的一道题，不强求满分做法，就看状压那部分吧。最关键就是想到有哪些因数能够影响我这一步的决策。动态规划应对多阶段决策问题时一般都需要这样考虑。仔细分析就会发现只有最近的(M)个才会有影响，然而题目又出乎意料的给了(M leq 5)，状压就很明显了
Code
注意矩阵存的是真正的状态，而不是状态的编号……
/*By DennyQi 2018*/
#include <cstdio>
#include <queue>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
#define int ll
const int MAXN = 100010;
const int MOD = 1e9+7;
const int INF = 1061109567;
inline int Max(const int a, const int b){ return (a > b) ? a : b; }
inline int Min(const int a, const int b){ return (a < b) ? a : b; }
inline int read(){
int x = 0; int w = 1; register char c = getchar();
for(; c ^ '-' && (c < '0' || c > '9'); c = getchar());
if(c == '-') w = -1, c = getchar();
for(; c >= '0' && c <= '9'; c = getchar()) x = (x<<3) + (x<<1) + c - '0'; return x * w;
}
int N,M,K,Ans;
int status[40],cnt,sta,_j,_k,exist[40];
struct Matrix{
int a[40][40];
inline void clear(){
memset(a,0,sizeof a);
}
inline void set_unit(){
memset(a,0,sizeof a);
for(int i = 0; i <= 32; ++i){
a[i][i] = 1;
}
}
};
inline Matrix mul(const Matrix a, const Matrix b){
Matrix tmp,res; tmp.clear(),res.clear();
for(int i = 0; i <= 32; ++i){
for(int j = 0; j <= 32; ++j){
for(int k = 0; k <= 32; ++k){
tmp.a[i][j] = (tmp.a[i][j] + a.a[i][k] * b.a[k][j]) % MOD;
}
}
}
for(int i = 0; i <= 32; ++i){
for(int j = 0; j <= 32; ++j){
res.a[i][j] = tmp.a[i][j];
}
}
return res;
}
Matrix ans,a;
inline bool check(int s){
int res(0);
while(s){
if(s & 1) ++res;
s >>= 1;
}
return res <= K;
}
void dfs(int x, int s){
if(x == M){
if(check(s)){
status[++cnt] = s;
a.a[(s>>1)+(1<<(M-1))][s] = 1;
a.a[s>>1][s] = 1;
exist[s] = 1;
}
return;
}
dfs(x+1, s);
dfs(x+1, s+(1<<x));
}
inline void ksm(int k){
while(k > 0){
if(k & 1) ans = mul(ans, a);
a = mul(a, a);
k /= 2;
}
}
signed main(){
N = read(), M = read(), K = read();
dfs(0, 0);
ans.set_unit();
ksm(N);
for(int i = 0; i <= 32; ++i){
if(exist[i]){
Ans = (Ans + ans.a[i][i]) % MOD;
}
}
printf("%lld", Ans);
return 0;
}